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Abstract. Co- limitation of marine phytoplankton growth by light and nutrient, both 
of which are essential for phytoplankton, leads to complex dynamic behavior and a wide 
array of coherent patterns. The building blocks of this array can be considered to be 
deep chlorophyll maxima, or DCMs, which are structures localized in a finite depth 
interior to the water column. From an ecological point of view, DCMs are evocative 
of a balance between the inflow of light from the water surface and of nutrients from 
the sediment. From a (linear) bifurcational point of view, they appear through a 
transcritical bifurcation in which the trivial, no-plankton steady state is destabilized. 
This article is devoted to the analytic investigation of the weakly nonlinear dynamics 
of these DCM patterns, and it has two overarching themes. The first of these 
concerns the fate of the destabilizing stationary DCM mode beyond the center manifold 
regime. Exploiting the natural singularly perturbed nature of the model, we derive 
an explicit reduced model of asymptotically high dimension which fully captures these 
dynamics. Our subsequent and fully detailed study of this model — which involves a 
subtle asymptotic analysis necessarily transgressing the boundaries of a local center 
manifold reduction — establishes that a stable DCM pattern indeed appears from a 
transcritical bifurcation. However, we also deduce that asymptotically close to the 
original destabilization, the DCM looses its stability in a secondary bifurcation of 
Hopf type. This is in agreement with indications from numerical simulations available 
in the literature. Employing the same methods, we also identify a much larger DCM 
pattern. The development of the method underpinning this work — which, we expect, 
shall prove useful for a larger class of models — forms the second theme of this article. 
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1. Introduction 

Phytoplanktonic photosynthesis provides the major biological component of the 
transport mechanism carrying atmospheric carbon dioxide into the deep ocean. 
Concurrently, plankton forms the basis of the aquatic food chain. As a consequence, 
phytoplankton growth and decay plays a crucial role in understanding climate dynamics 
[TO] and forms an integral part of oceanographic research. Conversely, climate changes — 
such as global temperature variations — have a direct impact on the aquatic ecosystem 
and thus also on phytoplankton [U [22]: there is a subtle and under-explored interplay 
between the dynamics of phytoplankton concentrations and climate variability At 
the same time, phytoplankton concentrations exhibit surprisingly rich spatio-temporal 
dynamics. The character of those dynamics is determined in an intricate fashion by 
(changes in) the external conditions, see [15] and the references therein. The building 
blocks for the observed complex patterns are deep chlorophyll maxima (DCMs) or 
phytoplankton blooms, in which the phytoplankton concentration exhibits a maximum 
at a certain, well-defined depth of the basin. These patterns are the manifestation of a 
fundamental balance between the supply of light from the surface and of nutrients from 
the depths of the basin. For the simplest models, in which spatiotemporal fluctuations 
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in the nutrient concentration are omitted (eutrophic environment), it has been shown 
that there can only be a stationary global attractor [17] . In particular, if the trivial state 
(no phytoplankton) is unstable, then there can only be a stationary globally attracting 
phytoplankton bloom with its maximum either at the surface (a surface layer), at the 
bottom (a benthic layer, BL), or in between (a DCM) [HI [12} [131 E] • This is no longer the 
case in coupled phytoplankton-nutrient systems (oligotrophic environment), although 
DCMs do tend to appear in those systems, also, for certain parameter combinations 
[TJ [TTJ US} CHI EE]- The detailed numerical studies reported in [15], however, show 
that the appearance of a DCM only triggers a complex sequence of bifurcations: as 
parameters vary, a DCM may be time-periodic, undergo a sequence of period doubling 
bifurcations, and eventually behave chaotically. 

In this paper, we focus on the effect that varying environmental conditions, and 
in particular nutrient levels at the ocean bed, have on the dynamics generated by 
the one- dimensional model for phytoplankton (VF)-nutrient (N) interactions originally 
introduced in H5l 

-VW z + [^P(L,N)-l]W, 

Y- l ^iP(L,N)W. 1 ' } 

In this model, the vertical coordinate z measures the depth in a water column spanned 
by [0, zb], while W(z, t) and N(z, t) are the phytoplankton and nutrient concentrations, 
respectively, at depth z and time t. As in [151 122], the system is assumed to be in 
the turbulent mixing regime [9} [13], so that the diffusion coefficient D is identically the 
same for phytoplankton and nutrient. The phytoplankton is characterized by its sinking 
speed V, its (species-specific) loss rate I, its maximum specific production rate /i, and 
its yield Y on light and nutrient. The model is equipped with natural no-flux boundary 
conditions at the surface for both phytoplankton and nutrients; the bottom is a source 
of nutrients but impenetrable for phytoplankton, 

DW Z -VW\ Z=Q , ZB = 0, iV 2 | 2=0 = 0, and N\ Z=ZB = N B . (1.2) 

The constant nutrient concentration Nb will act as the primary bifurcation parameter 
in this work. The nonlinear expression P(L, N) models phytoplankton growth due to 
light and nutrient, 

r(^) = {L + L ™ N + NH y (") 

in which Lh and Nh are the half-saturation constants of light and nutrient, respectively. 
(See [25] for a short discussion on the nature and specificity of P(L,N).) The light 
intensity L at depth z and time t is determined by the total amount of planktonic and 
non-planktonic components in the column [0, 2], 

L(z, t) = Lj e-Kiv'-RfS w{s,t)ds_ ( L4 ) 

Hence, the system is non-local — a typical feature of most realistic phytoplankton models. 
The light intensity term introduces an extra three parameters: Lj, the intensity of the 
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incident light at the water surface; K^g, the light absorption coefficient due to non- 
planktonic, background components and hence a measure of turbidity; and R, the light 
absorption coefficient due to plankton (self- shading). The first two of these parameters, 
together with zb, D, Y, and Ng quantify the effect that the environment has on the 
planktonic population. It is by varying these parameters that we examine the effect of 
changing environmental conditions on plankton. 

It is shown in [25] that the system (II. ip has a natural singularly perturbed nature. 
This can be seen by rescaling time and space via r = fit and x = z/zb and the 
phytoplankton concentration W, nutrient concentration N, and light intensity L via 

u > + {x,r) = ^-W{z 1 t) 1 rK*,r) = I-^M and j(x,t) = 

Substitution into ( 11. II) then yields, 

w r = ew+ - 2^u+ + (p{u+, ri, x) - i) 
Vr = £ (Vxx + l~ l p{u + , r), x) u + ) , 
with boundary conditions, 

(w+-2v^7Iu; + )(0) = (w+-2 v / V^u; + )(l) = and r) x (0) = r)(l) = 0. (1.6) 

For realistic choices of the original parameters of ( II. ip . 
D 
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cf. PJ3E5]. Effectively, e 1 ^ characterizes the extent of the zone where DCMs appear 
relative to the depth of the ocean. In this paper, we follow [25] and treat the parameter 
e as an asymptotically small parameter, i.e., we assume that < e 1 so that ( 11.51) 
has, indeed, a singularly perturbed character. The nonlinearity p in (I1.5P is given by 

p(u + , n, x) = r-j- — — — . - - — --, (1.7) 

with rescaled light intensity 

j(cu + , x) = exp I — kx — r / cu + (s, r)ds ) . (1.8) 



o 



The remaining six rescaled parameters of (jl.5|) . 

V 2 I . L H N H RDYNb 

V =T TT» * = -» 3H = ~j-, VH = 1 r r , K = K bg z B , and r = ,(1.9) 

4fiD n Li N B lz B 



are all considered to be 0(1) with respect to e in the forthcoming analysis (cf. 

Our attempt to comprehend the mechanism underpinning the appearance of 
phytoplankton patterns, as well as the character of such patterns, begins with the 
determination of the spectral stability of the trivial steady state u + = (0, 0) T . At 
that state, and in terms of the original system ( ll.ip . there is no phytoplankton — 
W(z, t) = — and the nutrient concentration remains constant throughout the column — 
N(z,t) = Nb, the value at the bottom of the basin (ll.2|) . The system (11.51) may be 
written compactly in the form 



er ha + et- 1 p(u+,T l ,x)u>+ 
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where 

,+ 



u 



'I 

Here, the nonlinear operator T + is densely defined in L 2 (0, 1) x L 2 (0, 1). The associated 
spectral problem has been investigated in full asymptotic detail in [25], where we worked 
with the linearization of f ll.lOp around u + = (0, 0) T , 



in which 



e d xx — 2 Jew d x + f — £ 
et x f ed 

f(*)= ,, U KX and "=TX— G (°» 1 )- ^ L12 ) 
1 + j H e KX I + Vh 

The spectrum a{VT + ) = {^ n }n>o U {A n } n >o of the operator VT + consists of two 
distinct, real parts associated with the two diagonal blocks of T>T + , cf. (11. lip . Here, 
the eigenvalues v n = —e(n + 1/2) 2 ir 2 are negative, independent of all parameters, and 
associated with the lower block. These eigenvalues, together with the corresponding 
sinusoidal eigenfunctions (0,cos((n + l/2)7rx)) T , describe nutrient diffusion in the 
complete absence of phytoplankton. It follows that the spectral stability of the trivial 
state is governed solely by {A n } n >o, the set of eigenvalues associated with the upper 
block. In [25], we identified two different linear destabilization mechanisms. In the 
regime v < /(0) — /(l), corresponding to reduced oceanic diffusivity or increased 
turbidity (cf. (II. 9p and (j!.12p ). the planktonic component u>q of the eigenfunction 
Wq associated with the critical eigenvalue A has the character of a DCM: Uq is 
localized in an 0(e 1 / 4 ) region centered around a certain depth at which it attains its 
maximal value, see Figure [TJ This depth can be determined explicitly: to leading order, 
/(x*) = /(0) — v [25]. Hence, x* increases monotonically from = to — 1 as v 
increases from v = to the transitional value v = /(0)— /(l). In the complementary case 
v > /(0) — /(l), corresponding to increased oceanic diffusivity or decreased turbidity, 
the planktonic component of the critical eigenfunction destabilizing the trivial state has 
the character of a BL: that is, it increases monotonically with depth and essentially all 
phytoplankton is concentrated in an 0(e 1//2 ) region over the bottom, see again Figure UJ 
In this article, we focus exclusively on the regime in which DCMs may appear, i.e., 
we assume throughout the article that v</(0) — /(l). In that regime, we investigate 
the nature of the bifurcation associated with the destabilization mechanism of DCM 
type. We know from [25] that, in this case, 

A n = A* - £^3^/3 jAi+i | + 0(£ l/2 ); (L13) 



with 



and where 



A, = /(0)-*-v=— - — -£-v (1.14) 

1 + JH 

a = F'(0) = -f(0)= K " 3H with F(x) = /(0)-/(x). (1.15) 

(1 + JhY 
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Here, A n < is the n— th root of Ai, the Airy function of the first kind. The bifurcation 
occurs as Ao crosses zero, yielding the bifurcation diagram in the left panel of Figure [2j 
More specifically, we focus on the (weakly nonlinear) dynamics generated by (11.101) for 
parameter choices such that 



where p > is fixed and A is allowed to be at most logarithmically large with respect 
to e. Note that one can tune the appearance of a destabilization of DCM type (i.e., 
of the simplest phytoplankton pattern) by choosing appropriately the parameters in 
( ll.lOP ; also, that Ao depends on all parameters with the exception of r, the rescaled 
self-shading coefficient, see the definitions of / and ao in (11.12ft and (11. 15ft . We remark, 
further, that the parameter v depends on the diffusion coefficient D (cf. f ll.9j) ). the main 
parameter varied in [15] and the one that most strongly depends on varying external 
conditions such as global temperature [22]. Finally, Ao is an increasing function of our 
bifurcation parameter N B through its dependence on u, see fll.l3j) — fll. together with 
the definitions of v in ( 11.12ft and of t]h in (11.9ft . Based on this final observation, we will 
often treat A as our bifurcation parameter. 

The first step in analyzing the dynamics generated by a linear destabilization 
mechanism is to perform a center manifold analysis to determine the local character 
of the bifurcation associated with the destabilization (see, for instance, [US]). This is 
a well-established procedure. In the setting of ( 11.16ft . this amounts to assuming that 
Ao is (asymptotically) smaller than all other eigenvalues, and it corresponds to the case 
p > 1 and A = 0(1). In this regime, the remaining eigenvalues {v n } n >o U {A n } n >i are 
negative and asymptotically larger than Ao, so that the local flow near the trivial pattern 
(0, 0) T is determined by the flow on the one- dimensional center manifold. The tangent 
space of this manifold at the trivial steady state is spanned by the critical eigenfunction 
Wq associated with Ao- Hence, this flow can be determined by expanding u + as 
u + (x, t) = £ p-1 / 6 f2o(r) Wq (x) + R(x,t), with f2 being an unknown, time-dependent 
amplitude and the higher order remainder R encapsulating the component of u + along 
directions associated with the stable eigenvalues — the additional 1/6 in the exponent of 
e follows from the projection analysis by which the equation for Q is determined (see 
below and Section |3]). An ODE for the unknown amplitude Qq is obtained through a 
projection procedure which is straightforward but can nevertheless be highly technical, 
especially in a PDE setting. In the case at hand, this equation reads 



to leading order. Thus, the procedure reveals the existence of a nontrivial fixed point 
which is stabilized through a standard, co-dimension one transcritical bifurcation. This 
fixed point corresponds to an asymptotically small DCM pattern, the amplitude of which 
grows linearly with A , 





— Aq flo — aooo(0) £l , 




u + (x) ~ e 



p - 1/6 Q* u+{x), withfi* 
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In general, one cannot expect to be able to compute the coefficient aooo(0) explicitly. 
Here, we exploit the singularly perturbed nature of (jl.lOp and the localized character 
of the eigenfunction Wq to do exactly that; in particular, it follows from the analysis to 
be presented in this article that 



see Section [31 In addition to yielding an explicit, leading order formula for the amplitude 
of the emerging (stable) DCM, this first result also implies that this DCM is ecologically 
relevant since the planktonic component of the primary eigenfunction is positive, Uq > 0, 
and hence also u + > by (IPty - fOB)) . 

The main aim of this paper is to develop an analytic approach through which 
one can go beyond the direct, finite-dimensional center manifold reduction outlined 
above. The original ideas underlying this approach — namely, the method of weakly 
nonlinear stability analysis — qualify as classical [23]. However, this particular method 
does not always provide more insight than the rigorously established center manifold 
reduction method: for instance, it also reduces the flow to a one-dimensional ODE 
of the form (11.171) . The situation is strikingly different here, as we can exploit the 
singularly perturbed nature of (11.101) . in conjunction with the asymptotic information on 
the eigenfunctions of T>T + obtained in (25] , to study in full analytic detail the case Ao = 
0(e) — see Section [4] — and even extend our analysis to the regime A = 0(elog 2 e) — see 
Section I4~5l This way, we can analytically trace the fate of the bifurcating DCM pattern 
well into the regime where the pattern undergoes secondary and, possibly, even tertiary 
bifurcations. 

For clarity of presentation, we divide the rest of the material in this Introduction 
into two parts. The first one focuses on the bifurcations undergone by the DCM patterns 
and on the ecological interpretation of our findings, while the second one focuses on the 
specifics of the asymptotic method developed in this work. 

1.1. The bifurcations of the DCM patterns 

The outcome of our asymptotic analysis is summarized in the right panel of Figure [2j 
The localized DCM that bifurcates as Ao crosses zero is a stable attractor of the flow 
generated by (11.11) . for all p > 1 and A = 0(1) with respect to e, cf. ( 11.161) . 
As we remarked above, the amplitude Qq of this localized DCM, and thus also the 
biomass associated with it, grows linearly with A in that regime, cf. ( 1 1 . 1 8 [) — ( fl . 1 9 p . 
Quite remarkably, from the point of view of our weakly nonlinear stability analysis, 
continues growing linearly with A also beyond the region where the center manifold 
reduction is valid. In particular, ( I1.18j) - (ll.l9j) remain valid in the regime p = 1 and 
A = 0(1), see ( 14. 9p . The corresponding biomass turns out to be 



Oooo(O) = (1 - - x*) 



a; /3 /(0)exp(|A| 3/2 ) 



>0, 



(1.19) 




1/2 



1 



(1 + 3h) 



u + (x) dx = e 



(1 - v) v{l - X*) 



A, 



o — 



(1 -v) (1 -x,) (£ + v) 



(1.20) 
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Figure 1. Left panel: a DCM profile for the planktonic component of (|1.5p (jl.OI) . 
Essentially all plankton is concentrated in an 0{e 1 ^) region around a finite depth x*. 
Right panel: a BL profile for the planktonic component of (|1.5|) - (|1.6|) . Here, essentially 
all plankton is concentrated in an 0(e 1 ^ 2 ) region over the depth of the basin. 



to leading order. This second result establishes that, in the Ao = 0(e) regime, the DCM 
pattern grows with v and hence also with Ng, the primary parameter measuring nutrient 
availability in the water column (see ( 11. 9ft and ( 11.121) ). This fact certainly reinforces our 
ecological intuition. 

The stability properties of the DCM mode corresponding to Qq, on the other 
hand, undergo a drastic change in that same regime. Our rather involved stability 
analysis of this emergent nontrivial steady state reveals that it becomes unstable, in 
this same Ao = 0(e) regime already, as A continues to grow and through a standard 
Hopf bifurcation; this is our third result. The appearance of this secondary bifurcation 
can be determined explicitly by our methods and, as we demonstrate, its onset occurs 
for values of A which increase unboundedly as x* j. (equivalently, as v j 0). It is 
natural, then, to attempt an extension of our analysis into a region where A 3> 1. In 
that regime, we establish the existence of a second localized DCM-type pattern: the 
associated reduced system has two critical points. Using our methods, we trace this 
second localized structure back to 0(1) values of A and find that it corresponds to an 
0{e 1 / 2 ) biomass depending nonlinearly on Ao- This is our fourth result. The stability 
type of this pattern can be also determined explicitly, although we do not undertake 
this task in the present work. 

Hence, our analysis yields that the stationary, stable, localized DCM pattern 
emerging at the transcritical bifurcation through which the trivial state becomes 
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unstable only persists in an asymptotically small, 0(e) region in parameter space 
before it yields to an oscillatory pattern emerging through a Hopf bifurcation. This 
fact reinforces our mathematical intuition that the appearance of this stationary DCM 
is the first step in a cascade of bifurcations leading to the chaotic dynamics reported 
in [15] — see also our discussion in [25]. In light of this, our analytical findings seem 
to agree qualitatively with these numerical results. In the same vein, our findings 
here suggest that the chaotic dynamics can be traced back to the small amplitude 
patterns emerging from the destabilization of the trivial steady state. (Of course, one 
must always exercise caution in interpreting numerical observations from an asymptotic 
point of view, especially when these simulations concern an unsealed system as is the 
case here: the authors of [15] have simulated the original system (11.11) and not the scaled 
system ( 11.51) .) Additionally, the fact that the onset of the Hopf bifurcation for v 4- 
occurs in the regime Ao 3> 1 — where certain higher order terms in our analysis become 
leading order and hence the analysis must be necessarily extended — possibly explains 
the absence of oscillatory and chaotic dynamics for small values of v, see [251 Figure 3.3]. 

Naturally, the questions on the fate of the oscillatory pattern generated through 
the Hopf bifurcation and on the nature of the larger DCM pattern are intriguing. At 
present, this is the subject of ongoing research. We do not pursue these questions further 
in this article, apart from a short discussion in its concluding section. 

1.2. The asymptotic method 

Parallel to understanding the character and fate of the linear destabilization mechanism 
established in [25], this article has a second — and from a mathematical point of view 
at least equally important — theme. Here, we have developed a powerful approach by 
which we can study the weakly nonlinear dynamics generated by ( 11.51) in full asymptotic 
detail and far from the region covered by more standard techniques (such as the 
center manifold reduction method). Indeed, one cannot hope in general to extend the 
analysis beyond the one- dimensional center manifold reduction discussed above and into 
the regime where Ao is not asymptotically closer to zero than all other eigenvalues. 
In other words, the sole analytical insight into the dynamics of the flow near the 
destabilization that one can generically obtain is the confirmation that DCMs indeed 
appear through a transcritical bifurcation. Let us look into this last point in more 
detail and for our specific model fll.5j) — (jl.6p . For Ao = 0(s) — equivalently, for p = 1 in 
( I1.16P — one can no longer 'project away' the directions corresponding to the eigenvalues 
v n = —e (n+1/2) 2 7r 2 associated with the operator VT + . Indeed, these are 0(e) for 0(1) 
values of n, and hence of the same asymptotic magnitude as A . As a result, the center 
manifold reduction approach yields a leading order system in at least asymptotically 
many dimensions. In general, such a system cannot be studied analytically, and one has 
to abandon the idea of performing an asymptotically accurate analysis. 

The crucial ingredient in our approach is our ability to explicitly determine, to 
leading order, all relevant coefficients in the reduced, asymptotically high-dimensional 
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Figure 2. Left panel: the bifurcation diagram for the trivial steady state of (|1.5j) in the 
regime v < /(0) — /(l) = F(l). The trivial steady state is stable in the region Ao < 
and unstable in the region A > 0. Here, udcm = ^(1 + in) and ubl — ^(1 + e K jif). 
Right panel: the bifurcation diagram for the small-amplitude DCM reported in (|1.17|1 
and (|1.19[) . The origin marks the transcritical bifurcation through which the trivial 
steady state is destabilized and the small-amplitude DCM pattern emerges. The value 
Aq marks the (first) Hopf bifurcation where this small-amplitude DCM is destabilized 
and a time-periodic DCM pattern is generated. 

system that extends the leading order, one- dimensional center manifold reduction. All 
of these coefficients are defined in a relatively standard manner in terms of projections 
based on the linear spectral analysis, see (I2.2ip in Section |2j We report the outcome 
of this part of our work in (14. ip . These leading order formulas clearly reveal a certain 
structure in these coefficients, which in turn reflects on the system of ODEs for the 
Fourier modes. It is this structure that allows us to extend our stability and bifurcation 
analysis. The sometimes remarkably subtle and laborious analysis by which these 
coefficients are computed provides the foundation for the strength and success of our 
program. Therefore, this analysis is a central component of our approach and lies at 
the core of the forthcoming presentation, see especially Sections [3] and EHZl 

An understanding of the conditions under which similar structure may be expected 
to appear is apposite to deciphering the fundamental mechanisms underpinning the 
success of our method and to determining a more general setting where this method 
is applicable. Naturally, what enables us to compute these coefficients, and thus 
also determine how they are related, is the accurate asymptotic control over the 
eigenfunctions that we establish. It is neither clear, a priori, that the structure present 
in the reduced system is a necessary consequence of that control, nor how much of that 
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control is necessary to establish the presence of sufficient structure. These issues are 
the subject of current research undertaken by the authors. Below we offer a brief sketch 
of the ideas behind this work in progress, as it also encapsulates the essentials of the 
method developed in the present work. 

To avoid the computational complexities associated with the weakly nonlinear 
analysis, we consider a much simpler, autonomous, coupled, reaction-diffusion system, 



Here, U and V are defined in [0, 1] x M + and obey certain boundary conditions, 
e.g., of homogeneous Neumann or Dirichlet type. The nonlinearities F(U,V) and 
G(U, V) are assumed to be smooth and at least quadratic in U and V; finally, 
< e C 1 is an asymptotically small parameter. The spectral problem associated 
with the trivial state (U, V) = (0, 0) decouples into two scalar problems of harmonic 
oscillator type. It immediately follows that, for v below a certain critical value 
v* , this trivial state loses stability when \x crosses a threshold /i*. Moreover, the 
eigenvalues {A^} n >o associated with the [/—component (and hence also with \x) are 
0(1)— apart, while the eigenvalues {A^} n > associated with the V— component (and 
also with v) are 0(e)— apart. Both sets of eigenvalues are naturally paired with simple 
trigonometric eigenfunctions. A straightforward center manifold reduction suffices to 
determine the nature of the bifurcation as fi crosses /i* and in the regime fi — fi* <C s. 
This situation corresponds directly to our — technically more involved — center manifold 
problem (|1.17j) - (ll.l9p briefly discussed earlier. Note that, here, the leading order analog 
of the DCM pattern identified in that discussion is a sinusoidal function. 

As long as fi — jj* <C e, the modes associated with the eigenvalues {A^} n 
remain slaved to the critical A^-mode, exactly as in our phytoplankton-nutrient model. 
However, this is no longer the case when \x — \x* = 0(e); in that regime, asymptotically 
many A^-modes are nonlinearly triggered by that critical mode. Nevertheless, the 
remaining A^-modes stay slaved, so that one obtains a reduced system of asymptotically 
high dimension. Here also, the coefficients of the leading nonlinearities can all be 
expressed in terms of projections along the eigenmodes, albeit they correspond to much 
simpler integrals. This process should enable us to study the conditions under which 
one is able to infer relations between these coefficients similar to those reported in ( |4.ip . 
This, in turn, should lead to conditions under which the reduced system has sufficient 
structure to allow a secondary bifurcation analysis — and perhaps even the identification 
of a cascade of subsequent bifurcations — of the nontrivial state bifurcating at \i — \f . 
An additional benefit of working in a simple setting of this sort is its amenability to 
rigorous analysis, which is beyond the scope of this article. 

A natural question to ask at this point is whether the model problem f ll.21j) 
shares enough structure with ( II. 5p to enjoy similarly complex yet tractable dynamics. 
Note, in particular, the absence of nonlocal and non-autonomous terms from (11.211) . 
Mathematically speaking, we expect these aspects to be insignificant for the type of 



U xx +fiU +F(U,V;e), 
e(V xx +uV +G(U;V,e) 



))• 



(1.21) 
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dynamics that the model exhibits close to bifurcation. (The situation is very different 
from the ecological point of view, naturally.) In the setting of (11 .5p . the nonlocality 
only complicates our analysis and thus clouds our understanding of the secondary and 
subsequent bifurcations beyond the center manifold reduction. Indeed, one expects the 
self-shading effect that a small DCM pattern has on itself to be much smaller than the 
shading due to the water column above it. This is most evident in Sections[3]and[6l where 
self-shading (quantified by the parameter r) is finally shown to contribute higher order 
terms only. Similarly, the sole role of the non-autonomous features of (II -ip is seemingly 
to introduce two spectra, {u n } n and {A n } n , with different asymptotic properties. In our 
model problem ( 11.211) . this is achieved instead by choosing disparate diffusivities for the 
two model components. 

Finally, it should be noted that our work resembles, but is certainly not identical to, 
Lange's work in pJ|J [20]. Lange has devised a powerful asymptotic method applicable 
to problems with closely spaced branch points which allows one to track the evolution 
of solution branches well into the regime where center manifold reduction breaks down. 
In our work also, the spectrum is asymptotically closely spaced, as are also then 
the branch points. Nevertheless, the differences between our work and the work in 
[T9] are substantial. Most prominently, Lange essentially defines branch points as 
points in parameter space where the linearization around the steady state admits a 
zero eigenvalue, see the derivation of [191 (3-10)] in particular. In our work, instead, 
the secondary bifurcation is induced by the parameter-independent negative spectrum 
related to pure diffusion and occur before any eigenvalues other than Ao have bifurcated. 
As such, these branch points are not captured by Lange's method. In fact, this secondary 
bifurcation — and, we expect, part of the cascade toward chaotic dynamics — occurs in 
a region of parameter space which is asymptotically small compared to the magnitude 
of the next critical eigenvalue Ai. Viewed from this perspective, then, the existence 
of the rich dynamics reported here for the regime Ao = 0{e) acts as a paradigmatic 
manifestation of nonlinear interactions. The linearly stable modes manage to have a 
decisive impact on the dynamics solely through nonlinear couplings and although a 
strictly linear point of view dictates that these modes should be utterly irrelevant. 

2. Evolution of the Fourier coefficients 

Our aim in this section is to write the PDE system (ll.lOp as an infinite-dimensional 
system of nonlinear ODEs and subsequently reduce it by relaxing the fast stable 
directions. To achieve this, we need explicit formulas for the (point) spectrum cx(T>T + ), 
as well as for the corresponding eigenbasis and its dual. The spectrum and the eigenbasis 
have been determined in [25]; we summarize the relevant formulas in Section |2J] below. 
We then obtain the dual basis in Section |2~21 by solving the eigenproblem for the adjoint 
operator {T>T + )* . Finally, in Section 12.31 we derive the desired ODEs for the Fourier 
coefficients close to the bifurcation point. 
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2.1. The spectrum and the corresponding eigenbasis ofVT + 

For completeness, we let H u + and H v be the subspaces of L 2 (0, 1) associated with the 
boundary conditions (ll.fip . H u be associated with the boundary conditions 

(d x U-y/^u){0) = {d x U-y^u){l) = 0, (2.1) 

and we write %+ = % w + x and H = T-L u x %r?- Both product spaces can be equipped 
with the inner product 

<«W> = ( ( ^ J > ( "I J \ = ("i + (*)" 2 + (*) + ^(x)»fe(x))dx. 

Subsequently, we define the function ^(x) = exp(^ v / e x) and the operator £ : % — > W+ 
corresponding to an application of the Liouville transform, 

(Eu\ (u+\ + ( u) \ ( uj+/E \ x + 

cm = I = — u with inverse u—\ = = £ u .(2.2) 

\ V J \ V J \V J \ V J 

(It is straightforward to check that the boundary conditions (II. 6ft for u yield the 
boundary conditions (12. ip for u + .) Both £ and £~ 1 are self-adjoint and bounded and 

VT = e->VT+S=^/-*- V E l), (2.3) 

with T>T densely defined and having self-adjoint diagonal blocks. 

As mentioned in the Introduction, the eigenvalues v n associated with T>T + 
correspond to the pure diffusion problem for the nutrient in the absence of plankton. In 
particular, they are solutions to the eigenvalue problem 

ed xx C,n = v n (n, with d x ( n (0) = Cra(l) = 0, 

and may be calculated explicitly, 

v n = -eN n , where N n = (tt/2 + nir) 2 for n > 0. (2.4) 

The corresponding eigenfunctions have a zero u + — component, and they are 





Cn 



where Cn{%) = \/2cos(\fNnx). (2.5) 



These are normalized so that HC1II2 = 1- 

The eigenvalues \ n , on the other hand, correspond to the eigenvectors 



w' 



Vn 

Here, the functions w+ and r\ n are solutions to 

ed xx u+-2^d x u+ + (f(x)-£-X n )u+ = 0, 
{d x u+-2y/^u+){0) = {d x u+-2yffiu+){l) = 0, 
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cf. (II. lip , together with the self- adjoint, inhomogeneous, boundary value problem for 
the component rj n , 

ed xx r] n - \ n r) n = -et^fu*, where d x r] n (0) = r/ n (l) = 0. (2.6) 

Equivalently, they are solutions to the self-adjoint, Sturm-Liouville problem 

e d xx u n + (f(x) -£-v-\ n )u n = 0, ^ ^ 

cf. (I2.2p — (I2.3p . As already stated, in [25] we derived the asymptotic expressions 

A„ = A* - e 1/3 a 2 /3 \A n+1 \ + 0(e 1/2 ), with n > 0, 

cf. fOB]) . Here, A* = /(0) - £ - v, <t = F'(0) = -/'(0), and A n < is the n-th root 
of the Airy function Ai, cf. (I1.15p . A formula for the n— th eigenfunction w n can also be 
derived using the WKB method, cf. [25J. The corresponding eigenfunctions for T>T + 
are w+ = (u^,r]n) T , where w+ = E u n — cf. (12. 2p . As we will see in the next section, it 
is natural to impose the normalization condition ||w n || 2 = 1. 

2.2. The dual eigenbasis ofT>T + 

To carry out the weakly nonlinear stability analysis of the bifurcating DCM profile, we 
also need to obtain the dual eigenbasis {w+} n > U {v n } n >o uniquely determined by the 
conditions 

w+) = (v n , v m ) = 5 nm and (w+, v m ) = (v n , = 0, 
for all n, m > 0. In this section, we show that 




Q n j and v n = | 

Here, u~ = u n /E, where u n solves the eigenvalue problem (12. 7p and satisfies the 
normalization condition 1 1 1 1 2 = 1. Further, expressions for the functions {( n }n were 
reported in (12. 5ft . while the functions } n may be found by solving the inhomogeneous 
problem 

ed xx iJ;- + 2,/evd x iJ)- + (f(x)-£-v n )il>- = -e^ -1 /Cn, (2q) 

^"(o)=^-(i) = 0. 1 • ] 

Alternatively, ip~ = ip n /E, where ip n solves the self-adjoint inhomogeneous problem 
e d xx i/j n + (f(x) - £ - v - v n ) i/j n = -e£~ l fE ( n , 

To verify the above, we start from the observation that the dual basis may be 
obtained by solving the corresponding eigenvalue problem for (T>T + )*, the adjoint of 
the operator T>T + . To calculate (T>T + )*, we write v~ = £~ 1 v, recall (12. 3p . and note 
that 

(VT + u + ,v~) = (VT + Eu,ir) = (£VTu,v~) = {VTu,£v~) = {VTu,v} = (u,VT*v). 
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This implies, further, that 

(VT + u + ,v~) = (u,VT*v) = (r^+.DT*^) = (u+^^VTSv-), 

whence (T>T + )* = £~ l T>T*£. Here, u + satisfies the boundary conditions (11.61) . whereas 
the boundary conditions for v~ are determined from v~ = £~ l v and the boundary 
conditions (12. ip for v — in particular, 

d x ijj- (0)=d x ^(l) = and d x ((0) = C(l) = 0, where v =1 ^ I. (2.11) 



It is straightforward to show that 
VT* = 



ed xx + f-£-v et-ifE 
ed xx 



and, since also (VT + )* = £- l VT*S, 

{VT+y = f ed xx + 2y/evd. + f-£ ei-^f ] ^ 

In view of (I2.12p . the eigenvalue problem (VT + )*w^ = A n w)+ for w+ = (&tiVn) T 
reads 



e d xx uj+ + 2y/evd x u+ + (f - £ - A„) w+ = -e£ 1 ffj n , 

subject to the boundary conditions (12. lip . The latter equation yields immediately 
fj n = 0, so that the former equation becomes homogeneous. It is now trivial to check 
that Cj^ = uj~ = u n /E, where u n solves the eigenvalue problem (12 .7p . This establishes 
the first part of 

Similarly, (I2.12p shows that the eigenvalue problem (VT + )*v n = v n v n has solutions 

Cn 

where the functions {^} n satisfy the boundary value problem (12. 9p . An application of 
the Liouville transform ip n = Eijj~ leads directly to the self-adjoint problem (I2.10p . 

2. 3. Evolution of the Fourier coefficients 

Our aim in this section is to write the PDE system (11. 1 [) as an infinite-dimensional 
system of nonlinear ODEs. We start by expanding the solution of d T u + = T + (u + ) in 
terms of the eigenbasis associated with the linear stability problem, 

u + (x, t) = £ C -V6 5 «»(t) w+(x) + e c J2 *n(r) v n (x), (2.13) 

n>0 n>0 

where c > is yet undetermined and the coefficients Q n and \& n are determined by 

n n = e- e 6- 1 {u + ,w+) and tf n = e~ c (u + , v n ). (2.14) 
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The exponent of 1/6 in the first sum of (12.131) is related to the localized nature of 
0Jq, the planktonic component of w^. In particular, Uq is shaped as a DCM with an 
O^ 1 / 6 ) biomass Hu^H (recall from our discussion following (12. 8 j) that, in contrast, 
| |cco 1 1 2 — 1)- More details on this issue will be presented in Section 14.3.21 Moreover, we 
have introduced the exponentially small parameter 

J — ( 



5 = exp 



<1, 



(2.15) 



the role of which is to counterbalance the exponentially large amplitudes of the 
eigenfunctions w£ and v n . In particular, 



J±(x) = y/v x ± I (x) and I(x) 



y/F(s) - F(x ) ds. 



(2.16) 



•I'd 



Here, the 0(e 1 ^ 3 )— parameter x corresponds to the turning point of (12. 7p 

xo = F-\\* - A ) = £ V V /3 | Ax | + O^ 1 / 2 ), 
while x* is the location of the DCM, the unique point where J_ 



maximum 



(2.17) 

attains its (positive) 

see also Appendix A ), i.e., 

= F-\v + F(x )) = F~\v) + 0{e 1/s ). (2.18) 

Thus, 5 _1 is a measure for the amplitude of the cu-component of the (linear) mode 
associated with a bifurcating DCM. The introduction of § in the decomposition (12.131) 
allows us to identify small patterns (u + <C 1) and is motivated by the observation that 
this decomposition yields 



(2.19) 



+ 

n ' 



V (x,t) = e c ~ 1/6 5E n >o^W^W + ^E n > *,WCnW 
The principal part of Uq is derived in Appendix A , while asymptotic formulas for uj, 
with n > 1, can be derived in a similar manner. For 0(1) values of n, it follows that cj+ 
is exponentially small everywhere apart from an asymptotically small neighborhood of 
where it attains its maximum value of asymptotic magnitude at most 0(£ -1 / 12 5 _1 ). 



Similarly, the principal part of rjo is given in |Appendix B[ together with an L°°— estimate 
which shows that rjo is at most 0(e 1 ^ 6 5~ 1 ) in [0, 1]. As a result, the coefficients of the 
eigenmodes f2 n (n > 0) in (12.131) are bounded uniformly in L°°(0,1) by an 0(e c_1 / 12 ) 
constant, while those of \& n (n > 0) are 0(1). In what follows, we derive the ODEs 
governing the evolution of these eigenmodes. 



2.3.1. Eigenbasis decomposition ofT + (u + ) To derive the ODEs for the eigenmodes, 
we need to express T + {u + ) in the eigenbasis {w^} n > U {v n } n >Q. In particular, we show 
that 



T + (u 



fc>0 



m>0 n>0 



W, 



k>0 



rank 



m>0 n>0 



(2.20) 
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where we have omitted an 0(e 3c_1 / 2 ) remainder. The coefficients appearing in this 
equation are given by the formulas 

a mn k = £~ 1/6 ( ( \a m u+,w£) = e~ 1/6 {a m u n ,u k ), 



1 

ei- 1 



a 'mnk = 8 1/3 ( „-i Sa m u+,v k ) = s 1/3 £ (a m u n , tp k ) + ei 1 (a m uj+,(k) 



1 



(2.21) 



b m u+,w^) = (b m u n ,u k ), 



ei' 1 

b' mnk = £ ~ V6 (^ J-i ) $t> m cj+,v k \ = £" 1/6 5[(6 m u; n ,^)+£r 1 (6 m a;+,a) 
Here, we have defined the functions 

a m = 5 [(1 - u)rj m + (1 - i/ -1 /)rs m ] /, with s n (x) = £ w+(s)ds, 

b m = {l-v)ftm. 1 ] 

Note that we use (•, •) to denote all inner products — in "H, and — as there is no 

danger of confusion. 

We start by decomposing T + (u + ) into linear and nonlinear terms by means of 

r+(u + ) = VT + u + + Af(u + ), where Af(u + ) = ^ \ (p - /) u + . (2.23) 

Substitution of the decomposition (12.131) into the linear term yields the eigendecompo- 
sition of that linear term, 

k>0 k>0 

= e c ~ 1/6 5 J2 x k tt k wt + e c J2 Vk*kV kt (2.24) 

fc>0 k>0 

where we have also used that w+ and v n are eigenvectors of VT + (see Section |2~T|) . It 
remains to express the nonlinearity A^(u + ) with respect to that same eigenbasis. First, 
since p — f contains the nonlocal term J* oj + (s)ds, see (ll.7l) - (ll.8l) . we write (cf. (I2.19P ) 

S(x, t) := e~ c+l/e / cu + (s, r)ds = 5 fi «.( r ) s n(z), (2.25) 

^ n>0 

where s n was introduced in (12.221) . We subsequently obtain, by (II .7p and (11.121) . 
1-r] 1 



V 



v 1 — 7] 1 + jn exp(nx) exp(e c 1 / 6 rS) 
1-n 1 



l-i/?7 1 + (1 - ^V)(exp(e c - 1 /6 ri S) _ I)' 

Substituting from (I2.19P for u + and r\ into this formula and expanding asymptotically, 
we find further 

P(C0 + , ^X) = f- 6^ J2 -S C J2 brr^rn + O (e 2 ^ 3 ) , (2.26) 

m>0 m>0 



— e 
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with a m and b m as denned in (I2.22p . We remark for later use that this asymptotic 
expansion remains valid for o(e 1 ' 4 ~ c ) values of VL n (n > 0) and o(e~ c ) values of *ff n 
(n > 0) (see our discussion following (l2~T9l) ). Next, (l2~T9l) and (12~26|1 yield 

(p - /) tu + = -e 2c -^ 5 E a ™< - 5 2c " 1 / 6 5 £ £ W m fi n , 

m>0 n>0 m>0 n>0 

where we have again omitted an 0(e 3c-1 / 2 ) remainder. By virtue of (12.231) . then, 

m>0 n>0 \ / 

~ 1/6 5 E E f J-i ) ^ + (^ 3c " 1/2 ) • 

m>0 n>0 \ / 

We may now decompose the spatial components in these sums with respect to the 
eigenbasis, 

I Ei~ l ) 6amU n = Y1 ( £l/66 amnk W k + £l/3 a 'mnk V k) , 
\ / fc>0 

! J 5 b m w+ = E i 6 bmnkwt + e 1/6 b' mnk v k ) , 

/ fc>0 

where the coefficients a mnk , a' mnk , b mnk , and b' mnk are found by means of (]2.2ip . Using 
this decomposition, we finally write (omitting throughout an 0(e 3c ^ 1 ^ 2 ) term) 

Af( U + ) = - e 2c E E E ( £ " 1/65 amnk W k + a 'mnk Vk) OA 



m>0 n>0 k>0 



- £2c E E E ( £ ~ 1/66 h ^k wt + b' mnk v k ) v m n n 

m>0 n>0 k>0 

= ~ S 2C ~ 1/6 5 E E E ( a mnk^m^n + b mnk ^ m fi n ) W k 
m>0 n>0 fc>0 

- ^ E E E ( a 'mnk^n + V k . (2.27) 

m>0 n>0 k>0 

Combining (I2.24p and (I2.27p . then, we arrive at the desired result (12. 20 p . 

2.3.2. ODEs near the bifurcation point We are now in a position to derive the ODEs 
for the amplitudes {Q n }n>o and {^ n }n>o- Differentiating both members of (12 . 131) with 
respect to time, we find 

d T u+ = e °-V6 5 E w+ + e c E^ ^ ( 2 - 28 ) 

fc>0 fc>0 

where the overdot denotes differentiation with respect to r. Next, d T u + = T + (u + ) and 
hence, combining (12.201) with (12.28)) . we obtain the ODEs for the amplitudes, 

(} k = \ k n k — e c E] E/ ( amnk ^"rrS^n + b mn k \I/ m f2 n ) + O (e 2c ) , 

m>0 n>0 

^k = ^k - e c E E ( a 'mnk + b' mnk V m Q n ) + O (e 2c ) . 

m>0 n>0 
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We now tune the bifurcation parameter A* so that the largest eigenvalue, Ao, is the 
only positive eigenvalue while the eigenvalues Ai, A 2 , . . . are negative. In particular, we 
write (cf. IOE]) ) 

A = eA , where < A << £~ 2/3 , 

v k = — eN k , where N k > is 0(1) for k — 0, 1, ... , 

A fc = - e 1/3 A k , where A k > for k = 1, 2, . . . . 

As we will see shortly, the cases of particular interest will turn out to be those where A 
is either 0(1) or logarithmically large. Note also that, since the distance between Ao and 
Afc is 0(e 1 ^ 3 ) by fll . 13|) . it follows that Ai, A 2 , . . . <C V\- Then, the evolution equations 
for the amplitudes become 

^mnO^m^n! (2.29) 

m>0n>0 m>0n>0 

V k = -eN k V k -e c J2H <nk^A ~ e c J2 E b 'mnk^A, k > 0, (2.30) 

m>0 n>0 m>0 n>0 

(l k = -e 1/3 A k Vl k - E C E E a mnkttA ~e C J2Yl h mnk^>A, k > 1, (2.31) 

m>0 n>0 m>0 n>0 

where we have omitted all higher order terms. 



3. Application of Laplace's method on a oo 

Explicit asymptotic expressions for the coefficients in the ODEs (12.291) — ( 12.311) obtained 
in the previous section can be derived by applying Laplace's method and the principle 
of stationary phase to the integrals in (12. 2 1 p . In this section, we demonstrate the use 
of the former in deriving an asymptotic formula for aooo- Asymptotic expressions for 
the remaining coefficients will be derived independently in Sections [5H3 after we have 
thoroughly analyzed the bifurcations that our system undergoes. Although the analysis 
in those sections is substantially more involved, our approach there is very similar to 
that in the present section. 

The main result of this section is the leading order approximation 

a 000 = A(A ) = aa(A ), (3.1) 

where we have defined the 0(1), positive, Ao— independent constant a and the function 
a by means of 

a = (l-u)f(0)C 1 C 2 a 1 /3 a;^>0 and a(A ) = ^j^ 1 Z^D . (3 . 2 ) 

V A cosh v A 

Here, a is defined in (1 1 . 1 5 [) . while 

/ roo \ -1/2 

Ci=yj Ai 2 (s)dsj , C 2 = exp(|A 1 | 3/2 ), and <r» = F'{x*) = -f'(x m ), (3.3) 
see [25] and Appendix A We start by recalling that the coefficient aooo is given by 

aooo = £~ 1/6 / a Q (x)u}%(x)dx, (3.4) 
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cf. (ESQ), where 

a ( x ) = 5 [r(l - u~ l f(x))s (x) + (1 - u)ijo(x)] f(x). 

Employing (12.221) . (12.251) . using the explicit approximation (IB. 51) for t]q from Appendix B 
and defining the functions 



we find further 

a (x) = 

Thus, 



hi{x,y) = f(x) 
fl- 



/(*) 



1 



V 



^0 



sinh 



%(x-y)) f(y) 



f(x) cosh (v^Ao^) 



f(y) sinh(v^(l-y) 



Ao cosh vAo 

: ^5 l h 1 (x,y)u+(y)dy + e- 1/6 5 I h 2 (x,y) cu+(y) dy. 



(3.5) 
(3.6) 



aooo 



i 

JO 



1 r l 



^0 



h 1 (x,y)ul(x)u£(y)dydx + e 1/6 5 
where Ii and X 2 are the two double integrals appearing in this expression. 



h 2 (x,y) ul(x)u^(y) dydx 



(3.7) 



We can obtain the principal parts of X\ and X 2 using Theorem Appendix D.2 
based on [23], in Appendix D We start with the latter integral which, as we will see, 



fully determines the leading order behavior of aooo- First, the normalization condition 
ll^olh = 1 yields h 2 (x, y) 0Jq(x) dx = h 2 (0,y) to leading order. Since, also, Uq has 
a unique maximum at the interior critical point x*, Theorem Appendix D.2[ I (with 
A = e" 1 / 2 , IT = - J_, and S = h 2 (0, ■)) yields 



Mo, j/) wo" (y) d y 



27r h 2 (0,x*] 



(3.1 



(^ 1/2 ) 1/2 

to leading order, where we have used the explicit leading order approximation ( 1A.2I) of 
cUq from [25] (see also Appendix A ), recalled the definition (12.151) of 5, defined 



and employed the identity J" = — 2~ 1 F~ 1 / 2 F' . 

Next, we show X\ to be exponentially smaller than X 2 . First, we rewrite it as 



(3.9) 



£ _l/4 C 2 Q-q y-v^ 



8vr 3 / 2 



hi(x,y) 



exp 



Hj(x,|/) 



(3.10) 



where we have used (1A.2I) and (A. II) . Here, D = {(x,y)\0 <y<x,0<x<l} and 



ni 




= J-(v) 


- 21 {x) 




and 


0i 


= 1, 


n 2 


'x,y) 


= My) 


-2/(1) 




and 


e 2 


= 29, 


n 3 


%y) 


= J-{y) 


+ 2I(x) 


-4/(1) 


and 


h 


= 6\ 


n 4 


[x,y) 


= My) 


- 2I(x) 


-2/(1) 


and 




= o, 


n 5 




= My) 


-4/(1) 




and 


h 


= 26 2 


n 6 


%y) 


= My) 


+ 2I(x) 


-6/(1) 


and 


06 


= e\ 



(3.11) 
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where I(x) and J±i(y) have been defined in (12.161) . and 



9 = ^— - ^= with o~ i = F(l). (3.12) 



Theorem Appendix D.l| yields, for each integral, a result proportional to 



exp(maxr Xj y) e DTlj(x,y)/y/E). We first identify max!!! and then show that maxl^ > 
max IT,-, for j = 2, . . . , 6; it follows that the dominant term in (I3.10P corresponds to Hi 
and the rest are exponentially smaller than it. Now, Hi has no critical points in D, and 
thus its global maximum lies on 

3 

dD = \J(dD)i = {(l,y)|0 < y < 1} U {(x,x)\0 < x < 1} U {(x,0)|0 < x < 1}. 

i=l 

First, the global maximum cannot be on (dD)i\ indeed, D lies to the left of (dD)i and 
d x Ili(x, y) = —2>JF(x) < 0, where we have introduced F(x) = F(x) — F(x ), so that 
Hi assumes higher values in D than on (dD)i. Next, Il 1 (x,x) = y/vx — 3I(x) on (dD) 2 , 
and thus maxll^x, x) = Hi(x**,x**) with < x** = F _1 (v/9) < a:* (recall (12. 18|) and 
note that F > is increasing). Finally, lli(x,0) = —2I(x) < on (dD)z, and thus 
max( 9D ) 3 IIi < < Hi(x**,x**). In total, then, we find that maxlli = Hi(x**,x**) > 0. 
Next, Il 2 (x, y) < Tli(x, y) < Hi(x** , x**) . Since the leftmost equality holds only in 
an 0(e ' )-neighborhood of x — 1, we find that maxn 2 < Tli( desired. 
Additionally, n 3 < n 2 on D, and thus also max^ILi < max^ri!. Next, n 4 has 
no critical points in D, and hence we need to examine its behavior on dD. First, 
the maximum cannot be on {dD)i by the same argument we used for II4. Next, 
n 4 (x, x) = J_(x) — 2/(1) on (dD) 2 , and thus max(aD) 2 n 4 = n 4 (x ! | t , x*) = J_(x*) — 2/(1). 
Finally, n 4 < -2/(1) < IT 4 ( x^ . x^ ^ on (dD) 3 , and hence maxll 4 = J_(x*) — 2/(1) = 
maxll2 < max 111, as desired. Finally, II5 < Il 4 and Ilg < II 4 , and the desired result 
now follows. 

These estimates show, then, that maxll! = H(x**,x**) > maxIL,-, for j = 
2,..., 6. Since (x**,x**) G dD and its Jacobian satisfies DHi(x**,x**) ^ 0, 
Theorem Appendix D.l yields for (13.101) the asymptotic formula 



Xl = c m c[ ^ CJCl exp (El^Zfl^ = e m cl exp {^fl) , 

for some 0(1) constants C[,C'-[ > 0. Since X 2 = 0(e l ' % 5~ 1 ) (EU) and, by (l2~T5|) . 
Xi _ 1/3 C? ( Hi{x*\x**)- J„{x* 

x 2 ~ £ c 3 exp { ^ 

with 

Hi(x**,x**) - J_(z») = [J_(x**) - J_(x»)] - 2I(x**) < 

(recall that x* is defined in (I2.18P as the location of the maximum of J_), it indeed 
follows that Xi is exponentially small compared to X 2 . 

We conclude that a oo is given by 5X 2 at leading order. Combining the expressions 
(13 .81) — ( I3T9T) with the definition of h 2 in (I3.6p . we obtain the leading order result (13. ip 
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by using the fact that /(a:*) = £, also at leading order. To derive this last identity, 
observe that — in the regime A <C 1 — it holds that A* = at 0(1), see fll . 14j) . f ll . 16j) . or 
equivalently that v = /(O) — £; further, and also to leading order, = v by ( 12. 18ft . 

so that the desired identity follows from the definition F(x) = /(O) — f(x) applied at 
x = x*. Finally, we note that higher order terms in formula (13. X h may be obtained solely 
by considering X 2 , as T\ is exponentially smaller than I 2 - 



4. Emergence of a stable DCM 

The trivial (zero) state is, by construction, a fixed point of the ODEs (1 2 . 2 9 h — (12.31 1) for 
the Fourier coefficients. In this and the next section, we identify the remaining fixed 
points of that system and determine their stability. In this entire section, we work 
exclusively in the regime p = 1 and A = 0(1). 

4-1. Asymptotic expressions for b m oo, a' 00k , and b' m0k 

As stated in the previous section, where we derived an asymptotic expression for a 000 , 
asymptotic expressions for the coefficients b m oo, Aqq^. , and b' m0k appearing in ( 12.29ft — (12.31ft 
are derived independently in Sections [5H7] below. Here, we summarize the leading order 
behavior of these coefficients, including also ( 13.1ft for completeness: 

a 00 o = A{Aq), 

bmoo = B, for m < e~ 1/3 , 

a'oofe = -^fc( A oM(A ), forO^A;<£ 

b'mok = -A' k (A )B, forO^A;,m«e- 1 /3. 

The function A was introduced in f[3TT]) - (l3T2]) . whereas B = V2 (1 - v) /(O) is a positive 
0(1) constant. Further, we have introduced the function A' k via 

A' k (A ) = a'a' k (A ), where ex = and a' k (A ) = ^^A _ u 2 ) 

Here, C 3 = (Ai'(Ai)) 2 . Note that, similarly to a (cf. 1ET2]I ). a 1 is an 0(1) constant 
independent of A ; the constants cr ; cr*, C\ and C 2 have been defined in ( 11.15ft and ( 13 .3ft . 
We also note the following identity concerning Airy functions (see Section 9.11(iv), 
identity (9.11.5)]) 

POO 

/ Ai 2 (s) ds = (Ai'iA^) 2 , or equivalently C\ C 3 = 1, 

which, in turn, yields an identity that will prove to be of exceeding importance in the 
rest of this section — namely, 

2a = a'B. (4.3) 

Asymptotic formulas for b m00 , a' 00k , and b' m0k and for higher values of m and k can be 
derived similarly. However, seeing as such formulas only contribute higher order terms 
in our analysis below, we refrain from presenting the details. In what follows, instead, 
we treat ( 14.1ft as being valid for all values of k and m. 
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4-2. The reduced system 

The system ( 12.29l) -( |2.3ip exhibits asymptotically disparate timescales depending on the 
value of p and associated with the asymptotic magnitudes of the eigenvalues. In this 
section, we investigate the case p = 1, in which regime Qq and ^o, . . . evolve on a 
slow timescale and the higher-order modes fii, Q2, ■ ■ ■ become slaved to them. Setting, 
then, p = 1 and rescaling time (with a slight abuse of notation) as t = er, the evolution 
equations become 

b mnQ ^ m n n , (4.4) 

m>0 n>0 m>0 n>0 

vj> fc = -JV fc * fc - e c ~ x E a ^ fi A - £ c - 1 E E Lk^, * > 0, (4-5) 

m>0n>0 m>0n>0 

e 2/3 n k = -A k n k - e c - 1/3 E a ™n k tt m n n - e c ~ 1/3 £ £ £w* m fi n , fc > 1. (4.6) 

m>0 n>0 m>0 n>0 

(Here also, the overdot denotes differentiation with respect to t.) It is natural to 
introduce slaving relations for the latter modes in this system, 

fi fe =s Ck G k (n , ^,#2,...), for all fc>l, (4.7) 

where the positive constants ci,C2, ... and the 0(1) functions (with 0(1) partial 
derivatives) G±, G2, ... are to be determined. To do so, we first write the evolution 
equations for Qq and ^i, ^2, ■ ■ ■ under these slaving relations; we find 

fio = Ao^o — £° 1 aooofio ~~ £ ° b m Qo^/ m , 

m>0 

9 k =- N k V k - £ C -V 00fc ^ _ £ c-i no b' m0k ^ k > 0, 

m>0 

where we have retained only the leading order terms from each sum. Dominant balance 
yields, then, c — 1. Next, the invariance equation for Q k yields that the right member 
of (14. 6p must vanish to leading order. Dominant balance yields c k = 2/3 and 



G k (Q , ¥ 2 , . . .) = -^nl - 5° v &m0fe * m 

At. ilfe 

K K m>0 



Recalling, also, (14. ip . we arrive at the evolution equations 

(4.8) 



tt = A Q - AQ 2 - BQ J2 m>o y m 



V k = -N k * k + A' k [An 2 + Bn E m > ^m] , k>0. 

Here also, we have retained only the leading order terms from each sum. 

Remark J^.l. The ODE (11.17j) — describing the flow on the one-dimensional center 
manifold in the regime where Ao = s p A <C e — can be obtained from the system (14.81) 
above as its A — > limit. Indeed, the \l/-modes become slaved to the mode Qq in this 
limit, and (14. 8p reduces to fll . lTj) with aooo(0) replacing A = aooo(Ao) (cf. (I4.ip ). Note 
that aooo has a removable singularity at zero, so we write aooo(0) = liniA ^o aooo(Ao) = 
(1 — £*) a. Using (13. 2p . it is plain to check that, indeed, the formula for a oo(0) reported 
in (II . 191) equals (1 — a;*) a. 
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4-3. The bifurcating steady state 

In this section, we identify the nontrivial fixed point of the reduced system 114. 8ft . In 
particular, we show that this fixed point is given to leading order by the formulas 

/ a \ A . TfVA x 2Ag cos( v / iVfca;*) 

where fc > and the parameter a was introduced in (I3.2p . Plainly, remains positive, 
and hence also ecologically relevant, for all positive values of Ao and all values of 
< x* < 1 (equivalently, all positive values of v up to the co-dimension two point). 
Further, the leading order expression 114. 9 j) for Qq exactly matches 

^S = — forA o^O, (4.10) 

cf. our discussion in the Introduction and in Remark 4.1 above. It will also be elucidated 
in Section 14.3.21 below that this fixed point corresponds to a DCM with an 0(e) biomass 
and an associated 0(e) nutrient depletion. 

Note that the denominators in the formulas for Qq and ^ k vanish for x* = 1. As 
explained in the Introduction, this value is attained by 2* at the co-dimension two point 
where DCMs and BLs bifurcate concurrently. This is another indication that the nature 
of the co-dimension two bifurcation is of independent analytical interest. 

4-3.1. Derivation of \4-9j ) First, setting the left members of 114. 8 j) to zero, we obtain 
an algebraic system for the nontrivial steady states, 

A -An -Bj2^m = 0, (4.11) 



m>0 



N k v k - A' k n 



0. (4.12) 



a n + b 

m>0 

Here, k > and we have removed a superfluous factor of f2 in 114. lip corresponding 
to the trivial steady state. Substituting from this equation into ( 14. 1211 . we obtain the 
equivalent formulation 

An + BY,^m = Ao and N k V k - A' k A fi = 0. (4.13) 

m>0 

This system is readily solved to yield 

fi„ = — and % = ^A n* (4.14) 

a'sBA + A k N k °' V ; 



where s is defined by the series 



A'm \ - 1 cos(viVma;* 

s 



a '^o N ™ ^0 + A o) 



To produce a closed formula for s, we recast this formula as 

E cosfy^i.) _L | cos (VNni x *) _ \ ^ cos( v / ]VmX ; .) \ . . 

n iV m (iV m + Ao) ~ A^ 2^ N m + A ' ' 

m>0 \m>0 m>0 
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with both series in the right member converging absolutely and uniformly with x*. The 
second series appearing in the right member of this last equation is a Mittag-Lefner 
expansion; analytic formulas for such expansions can often be obtained by means of the 
Fourier transform. In particular, [2TJ Eq. (1.63)] (with a = ir, b = iy/Ao, and I = 1) 
yields the explicit formula 

Ecos (-//Vm^*) _ sin (i\/Ao (1 - 2?*)) _ sinh (VA~o (1 - £*)) _ a(A ) , . 

iV m + Ao ~ 21^008(1^) ~ 2 v ^coshv / A^ ~ 2 ' 



m>0 

whence also 

Ecos (y/NmX^) _ a(0) _ 1 - 
iVm ~ ~2~~ ~ 2 ' ^ ' ' 

m>0 m 

Substituting into (j4.15p . we obtain 

1 - x* a(A ) 
S ~ ~2A 2A7' 
and therefore (14.141) for Qq becomes 

A 



(4.18) 



Q* = 

(a -a'B/2)a(A ) + (1 -x*)a'B/2 
The final formulas collected in ( 14. 9 p now follow by identity ( 14. 3 p and (I4.14p for ^>* k . 

4-3.2. Ecological interpretation We next proceed to show that the steady state 
(stationary pattern) we identified above corresponds to an 0(e) biomass with a 
corresponding 0(e) depletion of the nutrient. Indeed, ( 12. 19ft yields the leading order 
expression 

f tu + (x)dx = e 5/6 5tt* [ u+(x)dx (4.19) 
Jo Jo 

for the biomass. Here, we have also recalled that c = 1 and that f^, . . . are higher 

order, cf. (14. 7p . Recalling the definition of 5 in ( I2.15P and using the explicit leading 

order formula ( 1A.2jl for Uq, we obtain 

5 [ 1 toUx)dx = e-^^^ I" F- I l\x)^( j - {X) -J- {X * ) \ dx. 
Jo V 71 " Jo V V £ / 

As mentioned in Section El J-(-) has a sole, locally quadratic maximum at x*, and 

hence the integrand above is exponentially small except in an asymptotically small 



neighborhood of that point. Hence, the integral is of the type considered in Appendix D 



and Theorem |Appendix D.l| yields, to leading order 



o 2^ ^ F^(x*)^Jl(x* 

where we have also recalled that J"_ = —2~ 1 F~ 1 / 2 F'. Substituting back into (I4.19p . 
together with the formula for Qq given in ( 14. 9ft . we finally recover the first expression 
(ll.20p for the total biomass given in the Introduction. The second expression may be 
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derived by noting that (11.161) implies the leading order result vji\ + jjj) = £ + v, as well 
as that e Aq = u(l + ju)" 1 — £ — v. 

Similarly, ( 12.191) yields the leading order formula 

f r]{x)dx = e 5/6 5Q* [ r) {x) dx + £^2% [ ( k (x)dx. (4.20) 
Jo Jo fc > Jo 

Now, Ck(x) dx = (—l) k /Nk by ( 12.51) . Further, the integral J rjo(x) dx can be 
calculated using (IB. 1|) : integrating both members over [0,x] and using the boundary 
condition at zero, we find 



£A / ri (x) dx = £d x r} (l) 



f(x)u£(x) dx. 



(4.21) 



The derivative <9 x ?7o(l) can be estimated at leading order by (1B.5|) . Differentiating both 
members of that formula, we find 



tanh V A sinh V A (l -y)) - cosh yj A (l - y) f(y) ui+(y) dy. 



It follows from (14.2 ip . then, that 
£ Aq I r] (x) dx 



1 + tanh v 7 ^ sinh ^(1 - y) ) - cosh ^(1 - y) f(y) uj+(y) dy 



Applying Theorem Appendix D.l we obtain 



?7o(x) dx = e l ^5~ 



i CiC 2 tx 1/3 q~ 1/2 I cosh (ygj^sj 



An 



cosh vAn 



(4.22) 



which is the desired formula for rjo(x) dx. Recalling also (14. 9 p for ^/* k , we obtain from 
( I4.20p the leading order result 



where 
1 

s = — 

a 



n(x) dx = e ST q (Aq) 
A' m (A 



CiC 2 a o* ( cosh {t/Aq~x*) 



7 Ec- 1 )' 



N 2 



cos 



cosh v A 

(y/Nm.X*) 



+ a' s A c 



(4.23) 



iV£ (iV m + A ) 



v sin(y^:(l-x,)) 

A^(iV m + A ) " 1 " ] 



m>0 m m>0 m y u/ m>0 

This equation, together with ( 14. 9 p for Qq, yields the total nutrient depletion level to 
leading order. 



4-4- Stability of the small pattern 

In this section, we examine the stability of the DCM-like fixed point (f^,^*) = 
(Qq, \&o> • • •) which we identified in the previous section. In particular, we show 
that this fixed point is stabilized through a transcritical bifurcation at A = and that 
it subsequently undergoes a destabilizing Hopf bifurcation. 
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m>0 



* fc = -my 



A tt 2 + B tt 



m>0 



for all k > 0, 



around {Q* , Letting Q = Q* + dfi and # fc = + d^ fc and recalling (Q3j) . we 
find that the corresponding linearized problem reads 

dn = -An* dn - Bn* J2 d ^rn, (4.25) 

m>0 

d'^ fc = A' k [A + A ft*] dfi + [4 B fij - + A; 5 ft* Yl d ^™' ( 4 - 26 ) 

where we have only retained the leading order component from each term. 

Truncating at the arbitrary value k = K e N, we obtain the system <5$ = Lk 5$, 
where 5$ = (dft , d^ , d*i, • • • , d^^) T and 



/ 



K 



-An 







-5ft* 



-5ft* 



-BO, 



o 



\ 



A[ (A ft* + A ) A[ B ft* A[ BQq — N t 



A' BQ* 
A[ B ft* 



A' K (An* + A ) A' K B% A' K BQ* ...A' K Bn*-N K 

To characterize the spectrum of this matrix, we derive a formula for its characteristic 
polynomial det(£ — A/). First, we use the first row of £ — XI to eliminate the off- 
diagonal entries of all other rows. In this way, we find that the equation det(£ — XI) — 
is equivalent to setting to zero the determinant 



A + An* 



BVtt 



A'o(A-Ao) A + iVo 



5 ft* 




a;(a-a c 



A' K (\ - A c 



A + JVi 











5 ft* 




A 



Next, we can use the (k + 2)— nd column to eliminate the (k + 2)— nd entry of the first 
column, for < k < K, as long as A ^ —Nk- Since A = —Nk if and only if A' k = 
(as can be shown by expanding the determinant along the (k + 2)— nd row), we can 
eliminate all entries of the first column. (Note that A' k may indeed be zero: indeed, A' k 
is proportional to cos((fc + l/2)7rx*), which may or may not be zero depending on the 
values of k and x*.) Defining K, — {k : A' k ^ 0} C {0, . . . ,K}, KL k = fC — {k}, and 
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eliminating the entries of the first column as detailed above, we obtain 



Q(A) BQ* C 



A + N 
A + Ni 








X + N K 



0. 



(4.27) 



Here, 



Q(X) = (X + An*)l[(X + N k ) + Bn* (A -X)J2^k II ( X + N m)- 

fc£/C fcg/C m£K. k 

As detailed above, A = — Nk solves (I4.27P if and only if A' k = (equivalently, if and only 
if k ^ yVf). Further, A > cannot be an eigenvalue, since Q(A ) > and A + > 0, 
for all k G {0, . . . , K} — note that A, N , N±, . . . , Nk are all positive constants. Hence, we 
can extend the set over which we sum in the formula above to the entire set {0, . . . , K } 
and rewrite the equation for Q(X) in the form 

K 



Q(A) 



BQ* Q Y1 



Ai 



k=0 



N k + X 



X + Att* 

A - An 



(A -A)J](^ + A). 



As we just noted, the elements of the set {— N k } kG /c are not eigenvalues of C . 
Hence, the eigenvalues of £n are {— N k } k ^/c together with all solutions to 

K 



A'. 



k=0 



N k + X 



A + Att* 

A -An ' 



Substituting for A' k from (14. 2 p and for Qq from (14. 9p . recalling the identity (I4.3p . and 
letting K — > oo, we rewrite this equation in the form 

2A >^ cos( v / A 7 fcX*) A + A a(A )/(l - x*) 



E 



1 - ^ (N k + A) (N k + A ) 



A - An 



Here again, we may write 
cos (a/A^X*) 



E 

fc>0 



(N t + A) (N t + A„) A - A 



E 



cos 



fc>0 

1 a(A ) - a(A) 
2 



A fc + Ao 



E 

fc>0 



cos 



A fc + A 



A-A ' 

so that the eigenvalue problem becomes (1 — a;*) A 
a(0) = 1 — x*, we recast this equation as 



A a(A) = 0. Recalling that 



A 



a(0) 



-An, where we recall that a(A) 



sinh (Va(1 



. (4.28) 



a(A) a/AcosIia/A 
This equation is satisfied by some A if and only if it is also satisfied by its complex 
conjugate A*, as the right member is real and (A _1 a(A))* = (A*) _1 a(A*). Hence, we 
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may restrict arg(A) to lie in [0, n}. Further writing [i := a/A — fi R + i jjj, we rewrite the 
eigenvalue equation in its final form, 

. . (1 — x*)u 3 coshu 

pM:= - 3 mh((l-x.)ri A " «*«80.)elP.'/21- < 4 ' 29 ) 
We note here for later use that 

Re(pQLt)) , 2 2 . sinh[(2 - x*)//r] cos(x*/x/) - sinl^x*/^) cos[(2 - x*)//j] 
^RK^^i - V R ) 



x* cosh[2(l — — cos [2(1 — x*)^/ 

+ /i/(3//| - 



2 cosh[(2 - xA/i#] sin(x*/ij) - cosh(x*/i fi ) sin[(2 - x*)///] 



cosh[2(l - - cos [2(1 - x*)///] 

lm(p(fj)) 2 q 2 \ sinh[(2 - x*)fi R ] cos(x*/ij) - sinh(x*/z fl ) cos[(2 - x*)/x/)] 



x* cosh[2(l — x*)//r] — cos [2(1 — x*)//j 

+ /iii(3^-/i|) 



2 . cosh[(2 - x*)/x fl ] sin(x*/xj) - cosh(x*//,R) sin [(2 - x*)//j] 



cosh[2(l — x*)/i/?J — cos [2(1 — x*)///] 

^.^.S. Analysis of d4.29ty for A | We first establish that, as A | 0, the eigenvalues 
{A n }n>-i remain each in a neighborhood of the discrete values — Ao, — Ao, — Aq, . . .. 

For A = 0, (I4.29P yields either \i = (equivalently, A = 0) or cosh/i = (whence 
\i = iy/N m , m > or, equivalently, A G {— A" m } m > ). To investigate the possibility of 
negative eigenvalues A for A > 0, we set fi R ; = to find that (I4.29P reduces to 

p(i/i/) = r— ^ — ^ — r sin[(l - xA ///] cos /i 7 = A . (4.30) 

1 - cos [2(1 - x*)/i/J 

For A I 0, there is plainly a small root of this equation, /ii = \/Ao + 0(Ao), yielding the 
small eigenvalue A = — A + O(Aq). Additionally, all eigenvalues of the set {— N m } m >i 
perturb and remain real for small enough values of Ao > 0. Indeed, p(i-) intersects zero 
transversally at {y/N^} m >o, whence the persistence of any finite number of eigenvalues 
from among this set is automatically established. That the remaining, infinitely many 
eigenvalues also persist can be established by noting that, if the maximum value of p(i-) 
is positive between successive zeros, then this value grows unboundedly with fij. For 
the two first eigenvalues, in particular, we have the Taylor expansions 

A_i = -Ao + O(Ag) and A = -N Q + 4 ^1 Z x *) ^ Ao + O(A ), 

(1 - X*) 7T 

which demonstrate that both remain in the interval (—No, 0) and approach each other 
as Ao increases, see also Figure [31 These are precisely the two first eigenvalues that 
collide as A is increased, yielding a pair of complex conjugate eigenvalues. 

Next, the possibility of positive eigenvalues A — equivalently, positive solutions of 
(I4.29P — can be excluded by noticing that — A < while p(fi) > for all \x > 0. In fact, 
the possibility of eigenvalues anywhere but in a neighborhood of the negative axis can 
be similarly excluded by observing that 

1 /2 

i / \ i N i ,3 / cosh(2/x R ) + cos(2/x/) \ 

\ cosh 2(1 - x*)fjL R \ - cos 2(1 - x*)/z/ / 
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Figure 3. Plots of the function p(i/ij) (sec (|4.30[1 ) versus [ii (left panel) and versus 
A = — fj/j (right panel) for x* =0.7. Also plotted: the level line at p = Aq, here set 
at 0.2; the first two members of the sequence {\/Nk}k>o (left panel) and {— Nk}k>o 
(right panel) as solid dots; and the two smallest solutions fij to (|4.30[) (left panel) 
together with the first two eigenvalues A_i and Ao (right panel) they correspond to, 
all as hollow dots. 



Plainly, for every value of A , there exists a value fi* R (Ao) > which depends 
continuously on A , satisfies liniA -s.o /^(Ao) = 0, and is such that the equation 
= |A | cannot be satisfied for any hr > fi* R (A ). It follows that all solutions 
to ( I4.29P must lie in the half plane {fi\fif{ < fi* R (A )} which, in turn, corresponds to 
a neighborhood of the half axis {A G R | A < |/iJj(Ao)| 2 }. A local analysis around the 
origin now establishes the absence of eigenvalues with positive real parts, for A small 
enough, and hence also the result. 

4-4-3- Complexification of eigenvalues and the Hopf bifurcation As we briefly 
mentioned in the last section in conjunction with Figure El the two principal eigenvalues 
A_i and Ao come closer together as A increases. Eventually, they collide at a specific 
value fjfj G (0, 7r/2) and for A = A' = p(ifj?j) = max wg ( i7r / 2 ) pi}^i) > 0. For A > A' , 
this pair of eigenvalues becomes complex, so it is natural to examine whether it crosses 
into the right half-plane through the imaginary axis. (Note that no eigenvalues can 
cross through zero, as (14.281) does not admit a zero eigenvalue for A > 0.) 

To locate imaginary eigenvalues A = iA/ G iR, we set //r = /// = p, > and rewrite 
the real and imaginary parts of p as 

cosh[(2 — sin(:r*/i) — cosh(x*/2) sin[(2 — x*)jl] 
cosh[2(l — — cos[2(l — x*)jl] 
sinh[(2 — x*)p] cos(x*p>) — sinh^/i) cos[(2 — x*)p] 
cosh[2(l — x*)[i] — cos[2(l — x*)ji} 



Re(p(/x)) = 2(1 - x*)f 
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cosh[(2 — x*)p] sm(x*p) — cosh(x*/i) sin[(2 — x*)p] 



Im(p(/i)) = 2(1 - x*)p 6 



cosh[2(l — — cos [2(1 — x*)/j] 
sinh[(2 — x*)p] cos(x*p) — smh(x*p) cos[(2 — x*)p)] 



cosh[2(l — — cos [2(1 — x*)/j] 
The condition Im(p(/x)) = 0, derived from (I4.29p . yields 

cosh[(2 — x m )p\ sm.(x*p) — cosh(x*/i) sin[(2 — x*)p] 

= sinh[(2 — x*)p] cos(x*/i) — sinh^/i) cos[(2 — x*)p)]. (4-31) 

Therefore, the equation Re(p(/i)) = A , similarly derived from ( I4.29p . becomes 



3 cosh[(2 — x*)p] sin(x*p) — cosh^/i) sin[(2 — x*)fi] 



4(l-xJu d ^ *-tPi r 1 ! — =A . (4.32) 

v ^ cosh[2(l -x*)p] -cos[2(l -x*)p] V ; 

Condition (I4.3ip determines the values of p corresponding to imaginary eigenvalues 
A = 2i/Z 2 , while ( 14. 32ft yields the corresponding values of Ao for which these eigenvalues 
appear. Since the former of these can be recast as 

e (2-**)A sin (x^p - |) - e x ** sin ((2 - x*)p - ^) 

+ e- {2 ~ x *^ sin + J) ~ ^ sin (( 2 " z*)a* + ^) =0 ' ( 433 ) 

we see that there exists a whole, discrete sequence {pk}k>o of values p, see also Figure HJ 
As k — > oo, {pk}k>o limits to {(k + l/4)ir x^ 1 }^, the sequence of the set of zeroes of 
the first term in (I4.33P which becomes dominant in the regime p — > oo. Equation (I4.32p 
yields the leading order result 

A = 2v / 2ff 3 (l-x»)^" 3 (-l)Ve (w/4) ", as k -> oo, 

which establishes that the values pk corresponding to even values of k yield a positive, 
increasing sequence of values of A . (Odd k— values yield negative A — values.) In 
particular, the first Hopf bifurcation occurs at an 0(1) value of A when the complex 
conjugate pair (A_i,A ) crosses into the right half-plane through p . Higher, even 
A;— values correspond to Hopf bifurcations occurring at higher values of A , presumably 
when higher order eigenvalues cross into the right half-plane. 

These last remarks conclude our discussion of the DCM-like steady state for 0(1) 
values of Ao- In the next section, we investigate a logarithmic scaling for Ao in which 
the number of steady states of the system (I4.1ip - (l4.12p becomes two. 

4-5. A second DCM pattern 

So far, we have identified a DCM pattern corresponding to an 0(e) biomass which 
is stabilized through a transcritical bifurcation at A = and destabilized through a 
secondary, Hopf bifurcation at an 0(1) value of A . Here, we show that, the system 
( I4.4p - P~6j) admits a second, asymptotically larger, DCM-like steady state corresponding 
to an 0(e 1 / 2 ) biomass. We refrain from establishing the stability type, origins and 
eventual fate of that second steady state, reserving those problems to a later work. 
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Figure 4. Plots of the function in the left member of (|4.33[) for two distinct values 
of a:*, namely = 0.7 (left panel) and = 0.3 (right panel). The solid dots mark 
the members of the sequence {(k + l/4)7rx* -1 }fc>o. Note that the zeros of the plotted 
function approach this sequence rather quickly, with the quality of the approximation 
decreasing as i, f 1. Indeed, in that regime, all exponentials appearing in (|4.33[) 
remain approximately equal for a large range of p,— values, and hence the first term 
becomes dominant only in the far range. 



We start by noting that the inclusion of the first higher order term in the formula 
for a' 0Qk reported in (14.11) yields 

a' mk = -A' k (A )A(A )+e 1 / 2 aA' k (A ). 

This formula is derived in Section 1673] see ( 16. 2 j) in particular. Here, the A — independent 
constants a and a' were defined in (13. 21) and (14. 2ft . respectively, whereas the functions 
a and a' are reported in ( 13. ip and ( 14. 2p . Also, A' k (A ) = a! a'(A ) cos (v^V^x^), with 



1/3 -1/2 



a' 



a'(Ao) 



C1C2C0 ° 

V2f(0) ' 
sinh (a/A^(1 — 2*)) f^* f(x) cosh (y/A^x) dx 



(4.34) 



(4.35) 



V A cosh v A 

This formula for a'(Ao) is also valid in a logarithmic regime for Ao, see ( 1 6 . 3 f) for details. 
Since the first term in the formula for a' 00k above decreases exponentially with A (see 
( I3.2p ) whereas the second term decreases only algebraically, the two terms become 
asymptotically comparable for values of Ao logarithmically large in e, see Section [6] 
for details. 

Replacing the formula for a' 00k in ( 14. ip by the formula above, substituting into 
(I4.4p — (14.6]) . and working as in Section I4T2] we obtain the system 



fi 



A Q - 
-N k V k 



An 2 



(A' k A - e 1 * a A' k ) Q 2 + A' k B Q E m>0 



(4.36) 
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This is the analogue of (14. 8p with the inclusion of higher order terms. The fixed points 
(Qq, \&*) of this system are found, here again, by setting the left members to zero, 

A = Ao*+5E m > $; , (437) 



= -N k ** k + A' k A Q* - e 1/2 aA' k (n 



*\2 

o) ■ 



o- 



Solving the second equation for we find an explicit expression for in terms of Q 

Substituting this expression into the first equation in (I4.37p . we recover a singularly 
perturbed quadratic equation for 

£l/2aB f E 7?) W - f A + B A ^ ) n* + A = 0. (4.38) 

Vm>0 iV "V \ m>() "V 

In Section 14.3. 1[ we obtained the formula 
A + BA ^^ = aa(0), 

m>0 m 

while (I4.16P yields 

^^- = «a(A )^ — = _ aa(Ao). 

m>0 m>0 

It follows that the quadratic equation yielding Qq can be recast as 
1/2 &'Ba>(A ) JV^ _ 

£ 2 ™ iZo+ aa(0)" U - 

The two solutions of this equation are 



o*,± = - 1 / 2 l±v / l-2 g 1 / 2 a^Ao^(Ao)/( a q(0)) = ; 2e~ 1 / 2 / (a 7 £? a^Ap)), 
a'5S'(A ) I A /(aa(0)), 



with the first one corresponding to the asymptotically larger DCM-pattern and the 
second one corresponding to the DCM-pattern identified trough our earlier work. We 
remark here that this first steady-state is, indeed, within the reach of our asymptotic 
methods, as Qq and safely remain asymptotically smaller than the asymptotic bounds 
£ -3 / 4 and e" 1 for which our work in Section 12.3.11 remains valid. Note also that this 
steady state is a nonlinear function of A , with the distinguished limits 

2 e -i/2 4e" 1 /2 

lim fi*(A ) = - -g. and tt* (A ) = — — A , as A oo. 

A ->o (1 — x*)a'B J Q j(x)dx la' B 

In particular, this second pattern approaches a nonzero value as A 4 and eventually 
grows linearly for Aq 3> 1. 



Emergence of localized structures in a phytoplankton-nutrient model 
5. An asymptotic formula for 6 m0 o 



34 



In this section, we derive the asymptotic formula for b m oo given in (14.1jl . where m G N 
and 

«i 

Km = (1 - v) / f(x) ( m (x) Uq(x) dx. (5.1) 



o 



As detailed earlier, the function u , appearing in (15. II) . decays exponentially outside an 
0(e 1 ^ 3 ) — neighborhood of the origin (cf. ( 1A.1I) ). whereas the period of the sinusoidal 
term ( m is equal to 2ir/\ / N m = 4/ (2m +1). Below, we analyze the three different 
regimes — in which the integrand is predominantly localized, concurrently localized and 
oscillatory, and predominantly oscillatory — and we derive the leading order, uniform 
asymptotic expansion 

b, for m < e~ 1/3 , 

£l/3 Ai ' (? + ^) for m = °^" 1/3 )' (5.2) 



J m00 



bCf I 
'o 



Here, b = y/2 (1 - v) f(0), cf. Section gU 
5.1. The case m <C e" 1 / 3 

Here, 2ir/ y/N m 3> e 1 ^ 3 and hence the integrand is predominantly localized around x — 0. 
Thus, Cm m ay be approximated to leading order by £ m (0) = \/2 in that neighborhood. 
Since 1 1 Wo || 2 = 1 (cf- our discussion in Sections l2.lH2.2p . we obtain the desired formula 

b m00 ~ b. (5.3) 



5.2. The case m = 0(e 1 ^ 3 ) 

Here, 2n / 'y/N m = O^e 1 ^ 3 ), and hence the neighborhood of the origin outside which coq 
decays exponentially and the period of the sinusoidal term are of the same asymptotic 
magnitude. Defining the new variable r = T\ x in 1 15. ip . with T\ = \A\\ /x 1 12. 1 7[) . we 
obtain 



V2(l 



■>mO0 



/ I — J cos 



n Jo V r i/ V n 

Now, (lA.ip yields, to leading order and for any r <C £ _1 / 3 , 



(5.4) 



e" 1 / 6 d a 1/6 Ai (r + Ax 

2 v ^FF 1 /4( T /ri) eX P 



for r G [0, -Ai), 



■ n 



F(x ) di 
for r G (-Ai,r ], 

where we have also changed the integration variable by means of s = t/ri+x^. These two 
formulas agree — as, indeed, they should by construction — in the regime 1 < r < e^ 1 ^ 3 . 
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Indeed, recalling the asymptotic expansion of Ai in a neighborhood of infinity [2], we 
find that the first branch of the formula above yields 



-1/6 1 °0 



exp (-^(- + ^i) 3/2 )- 



2 v ^r 1 /4 

Similarly, the formula in the case r G (|Al| ,tq] becomes, upon Taylor-expanding F, 
E ' o /i.im ex P 



2^rV4 ~^ ^ 3 V e \ r x 

That the two formulas agree now follows from the definition T\ = \Ai\/xq and the 
formulas (12.171) and (13.31) for x and Ci- Hence, we may write 



C^j ~ ^ 1/6 a 1/6 Ai (r + A x ), for r<e 



1/3 



Since the contribution to the integral in (15 .4p of greater values of r may be estimated 
to be exponentially small, we can write 

b m00 = e^ 3 ^ l -») C l^ f (L) cos (V^) Ai 2 (r + A,) dr 



oo 



( el/3 ^F r ] Ai2 ( 







-1/3 



to leading order, as desired. Note that this formula reduces to (15. 3p . for m <C e 
5.3. The case m ^> e^ 1 ^ 3 

Here, 2n/ y/N m <C e 1 ^ 3 . Similarly to our work in the previous section, we define the new 
variable r = e~ l l 3 x. We find, then, 

b m00 = V2e 1/3 (1 - v) J g{r) cos (V /3 v^r) dr, 



where g(r) = f (e^r) loq (e^r). Using Theorem Appendix D.4 (with A = e 1 ^ 3 y/N, 



$(£) = t = r, and h(r) = g(r)) and the fact that the right-boundary term is 
exponentially smaller than the left one, as ojq(1) is exponentially smaller than ujq(0) 



(cf. Appendix A ), we obtain 

fc+i 



b m00 = V2e^ 3 (1 - u) Re (E^O) {^/^ ) 

-V2^E(- 1 rV-'(0)(^) . (5.6) 



Recalling the definition of g, and employing ( lA.lj) and that Ai(Ai) = Ai"(Ai) = 0, we 
calculate 

</(0) = and g"'(0) = -6[Ai'(A 1 )} 2 C 2 1 al 
The desired result now follows, while (15 .6p also reduces to (15. 5ft for m = 0(e^ 1 ^ 3 ). 
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6. An asymptotic formula for a' 00k 

In this section, we derive the asymptotic formula for a' 00k for 0(1) values of Ao collected 
in flED, 

a' 00k = -A' k (A ) A(A Q ), for ^ k « e^. (6.1) 
Further, we extend this result to 

«oo fe = " (4(Ao) ^(Ao) - s 112 a 4(A )) , for ± k « e" 1 / 3 , (6.2) 
which remains valid at least in the regime 

A = — -J log 2 £ + log e log(- log e) + —J log 2 log e + // log e, (6.3) 

J^i-t- ^ J- ^ >X>;f; 

for all ji G (— oo, jt/o] an d yUo > any 0(1) value. Here, A(A ) = aa(A ), A' k (A ) = 
a'a' k (A ), and A k (A ) = a'a'(A ). The A — independent constants a, a', and a' were 
defined in (13. 2p . (14. 2p . and (14 .34 p . respectively, whereas the functions a, a', and a' 
are reported in (13. ip . (14. 2p . and (I4.35p . We remark, here, that these results are only 
valid for those values of k for which Ck( x *) 7^ 0- F° r the remaining values of k, 
Theorem Appendix D.l yields (algebraically) higher order results. Also, we note that 
asymptotic formulas for higher values of k can be derived as in the previous section, 
albeit at considerable extra computational cost. 

We first write out explicitly the expression for a' 00k yielded by (I2.2ip . 

a' 00k = e~ 1 ^ 3 5 / a (x) uq(x) ip k (x) dx + e 2 ^ 3 5£~ 1 / a (x) (x) ( k (x) dx. 
Jo Jo 

Recalling the definition of ao from (I2.22p and working as in Section [3j we obtain further 



4 00fc 



e 



1 r x 



JO 



' 1/3 S 2 / / h 1 (x,y)u^(y)uJo(x)ip k (x)dydx 



+ e 1/3 5 2 / / h 2 {x,y)uj+(y)uj Q (x)%l) k (x)dydx 
Jo Jo 

+ e^i- 1 [ f h(x, y) ( k (x) u+{y) co+(x) dydx 



JO 
1 pi 



JO 



+ e 2/3 5 2 r 1 / h 2 (x,y)( k (x)u+(y)u+(x)dydx. 



Substituting, finally, from ( 1C.2ip . we obtain an integral formula for a' 00k which is 
amenable to the sort of asymptotic analysis employed in Sections [3] and [51 
e-l/302 r rl rx rx 



a 00k 



(W^) 1 h 1 , k (x,y,z)u (x)i/j k -(x)uj h (y)i> k h t+ (z)dzdydx 
Jo Jo Jo 

rl rl rx 

+ (W^,y 1 / / h 2 , k (x,y,z)uJo(x)ilj k -(x)uj£(y)tp~ k h + (z)dzdydx 
Jo Jo Jo 

rl rx rl 

- (W^)" 1 D k (ti) I \ \ h hk {x,y,z)oj {x)i) k ^(x)ujQ{y)il)l_(z)dzdydx 

Jo Jo Jo 

- (W^y 1 D k (0) / / / h 2 , k (x,y,z)u}o(x)ip k -(x)uj£(y)'ip£_(z) dzdydx 

Jo Jo Jo 
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Figure 5. The domains of integration for the integrals Zi, . . . ,Xg in (|6.4p . 



pi px pi 

+ ( w ^ 1 / / h 1>k {x,y, z)u {x)if; k>+ {x)u^{y)ilj^_{z)dzdydx 

JO JO Jx 

+ (Wip)' 1 / / / h 2 , k (x,y,z)u (x)i{; k>+ (x)u£(y)i(;^_(z) dzdydx 

Jo Jo Jx 

pi px 

+ e / h 1 (x,y)Ck(x)u)£(x)u)£(y)dydx 
Jo Jo 
1 pi 



JO 



+e / h 2 (x,y)Ck(x)u£(x)u£(y)dydx 



(6.4) 



Here, hi tk (x, y, z) = hi(x, y) f(z) Ck(z), for i — 1, 2, and the constants D k (0) are reported 
in (lC.19j) - (jC.20p . Let X l7 . . . ,X 8 denote the integrals in the right member of H6.4[) in the 
order that they appear (the three-dimensional domains of integration for T\ y . . .,Xq are 
sketched in Figure [5]). In what follows, we omit the term 9ojq_(1; xq) uq )+ (x; xq) in the 
expression flA.ip for uq, as one can show that its contribution is exponentially small 
compared to the leading order terms (see also Sections [3] and E}. 

6.1. A rewrite of flff.^p 

In this section, we group together integrals appearing in the right member of (16. 4p in 
order to achieve a first reduction in the numbers of terms of that member. We start 
with rewriting the term -(W^)" 1 _D fc (0)Z 4 + {W^y 1 X 6 + el 8 . First, 

X 4 = / / (/ h 2 ,k{x,y,z)u (x)$ k -(x)dx\w£(y)i(>£ (z)<LA ys , 

J Jtt x D 4 \Jo J 

where ir x is the orthogonal projection on the yz— plane — and hence ~k x Da = [0, l] 2 — and 
dA yz is the area element on that plane. Since i\) k - = e 1 ^ 3 Ci 1 a 1 Uq in a neighborhood 
of the origin (cf. ( lA.ip and ( 1C.12P ). ujq is exponentially small outside this neighborhood, 
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and 1 1 tcJ 1 1 2 = 1; we write 



X 4 = £ 1/3 Cr 1 a 1/3 / / h 2 , k (0,y,z)u+(y)i;+_(z)dA yz . 

Recalling that ip^ _ = Eip k _, according to our convention in Section [2, and substituting 
into the formula above from (1A.2|) and (1C.12I) . we obtain 

whence, employing also ( 1C.19I) . we find 

(W^ 1 D k (0)X 4 = £~ 1/6 J jf ^ E 4 (y,z) exp (^=^j dA yz . (6.5) 

Here, 

„ C\d k h 2 {0,y) f{z)( k (z) 

" 4(y '* ) = 4^ F^(y)F^(z) MV:z) = J-(y) + J-(z)- (6-6) 

Next, we rewrite X 6 , 

(W^y 1 l e = (W i ,)~ 1 I! ( I h 2 ,k(x,y,z)u (x)ilj k>+ (x)dx) u£ (y) ip£ t _(z) dA 

J Jtt x D 6 \Jo 

Employing flA.ip and (IC.13j) . now, we obtain 

m- 1 * = - 1/6 ^ / / ( r h2A ^ z) dx i tuz) dA yz . 



(6.7) 



yz- 



2ttW^ J J TxD6 \J y /W{ Xr 
Further using (1A.2I) and, once again, ( 1C.13I) . we find 



{W^X, = £ V3 J J H6(y? z) exp (5^) dA yz . 



Here, n x D 6 = ir x D 4 , U 6 (y, z) = U 4 (y, z), and 

Similarly, renaming (x,y) as (y, z) in X 8 , we derive the formula 

els = J jf E 8 (y, z) exp {^^) dA yz , (6.9) 

where D 8 = ir x D 6 = n x D 4 , U 8 (y, z) = U. 6 (y, z) = n 4 (y, z), and 

C\CloT h 2 (y,z)( k (y) 
^ Z) = ~^r~ FV*(y)FW{zY (6 ' 10) 
Combining f l6.5p - fl6.10p . we obtain 

- (wg- 1 D k {o) x 4 + (wg- 1 x 6 + x 8 

= -^ 1/6 J jf ^ S 4 (y, z) exp (^|^) ^ (6.11) 
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where, to leading order, uniformly over -k x D±, and for all 0(1) values of A 



o- 



~, x x Cf 4 h 2 (0,y)f(z)( k (z) 

My, z) = My, z) = 47r ^ FifiMFW® • (6 ' 12) 

Next, we rewrite the term {W^)^ l X 5 + eX 7 . We write first 



J Jtt x D 5 \Jy J 

where 7r XJ D 5 = {(y, z)|0 < y < z, < z < 1}. Now, flA~T]) - flA~2|) and f]CTP2f yield further 

= £ " 3 / jf ^ Ssfez) exp (^5^) <M,.. (6.13) 
where we have defined the functions 

c;cj g „ 2/3 /MOM f tiM , ,„,., 

^• 2) = 8^, FV*(y)FV\z) J, 7^ ( " 4) 
and n 5 (?/, z) = n 4 (|/, z). Next, renaming x into z in X 7 , we find 

el 7 = j jf 3r(y, z) exp {^^) dA yZ} (6.15) 

where 

n n \ CfCia 2 /3 hi(z,y)( k (z) . , , w . 

D 7 = tt x D 5 , &r(y,z) = — jri/4( y ) jri/4( z ) ' and n 7(l/^) = n 4 (y, z). (6.16) 

Combining f l6.13p -f l6.16p . we find, to leading order and uniformly over D 8 , 

(W 4! )- 1 l 5 + el 7 = e 1/3 j jf Z 5 (y,z)exp(^^j dA yz , (6.17) 

where E 5 (y, z) = E 5 (y, z) + e 1/2 E 7 (y, z). 
We now rewrite (W^)~ 1 Z 2 . First, 

1 2 = H 2 (x) f(z)( k (z)uo(x)ip k ^{x)ipl + {z)dA xz , 

J J -KyDl 

where H 2 {x) = h 2 (x,y) (jJq{v) dy. Substituting for oj^{y) from (I A . 2 [) . we find further 

X 2 = e~ l/l2 Cl ^ /3 [ [ H 2 (x)uj (x)^_(x)f(z)( k (z)i;+ + (z)dA xz , 

V 71 " J Jtt v D 2 

where H 2 (x) = L h 2 (x,y) F~ l / A {y) exp( J_(y)/y/e) dy. Using Theorem Appendix D.l 
now, we obtain 



(W^r 1 X 2 = e 1 / 6 r 1 C£'/ / h 2 (x,x*)u; (x)i; kt 4x)f(z)( k (z)i;+ + (z)dA xz 

J Jn y D2 

= e llX2 j jf d E 2 (x, z) exp (^^) dA xz , (6.18) 



Emergence of localized structures in a phytoplankton-nutrient model 40 

where C 2 is an 0(1) constant, 7r y D 2 = {(x,z)\0 < z < x, < x < 1}, U 2 (x,z) = 
J-{x.) + J+{z)-2I{x), 

E 2 (x, z) = C 2 M^AJ^l^A , with C 2 an 0(1) constant. (6.19) 
V-F(x) F 1 ' 4 (z) 

Finally, we rewrite (W^,) -1 D k (0) X 3 . First, 

2a = ( / f( z )Ck(z)ipk-(z)dz) / / hi(x,y)u (x)i/)k,-{x)u£(y)dA ai y. 



Substituting from (1A.1I) - (1A.2|) and (1C.12I) into this formula and interchanging the roles 



of y and z in the single and double integrals, we find 




U 3 (x,z) 



(W^)- 1 D k (0)X 3 = e-^~h / / %{x,z) exp -**pL dA xz , (6.20) 



TTyD2 



where I 3 = F 1 / 4 (y) f(y) ( k (y) exp (J_(y)/y/e) dy and 



E 3 (x,z) = C 3 - ^^] and IL 3 (x, z) = J.(z) - 2J(x), (6.21) 
for some 0(1) constant C3. 

j4n asymptotic estimate for a' QOk in the regime A = 0(1) 

In this section, we estimate the various terms derived above, starting from 
-(W,^ 1 D fc (0)Z 4 + (W^)~ 1 X 6 + el 8 (cf. flBUMBIEJ)). The exponent n 4 becomes 



maximum at the interior critical point (x*,x*), and thus Theorem Appendix D.l yields 



-(W^)' 1 D k (0)X 4 + (W^)" 1 X 6 +X 8 = -^ 1/2 ^T^ (V 1/6 H 4 (x*,**) exp H^^A Yj 
where 

C 4 = C\ (o-jW^)" 1 d k Ck(x*) /(a;*) /i 2 (0, 2*). 

Next, we estimate (W0) _1 X 5 + £:X7, cf. ( 16. 17(1 . The sole (quadratic) maximum of II4 
in D 7 lies at the critical point (x*,x*) G <9_D 7 , where H 5 (x*, x*) = and ^(x*, x*) 7^ 0. 



Recalling the definition of H5 and employing Theorem Appendix D.l then, we obtain 



(W^Xs+Xy = e (e^Co exp ( ^^T^ )) = ^r 2 ^, 

for some 0(1) constants Cq and C7. 

We now estimate the remaining three integrals starting with (W^)~ 1 X 2 , cf. ( 16 . 1 Sj) — 
( 16 . 1 9j) . The exponent U 2 has a sole maximum at the point (x*,x*) G d(n y D 2 ) which 
is not a critical point (compare to the maximization of II4 in Section [3]). Hence, 



Theorem Appendix D.l yields 



(W^X 2 = e^C 2 " (V 2 ~ 2 (x*,x*) exp ( ^^f^ )) = £ 4/3 <T 2 C 2 
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for some 0(1) constants C' 2 n and C 2 . Next, since D\ C D 2 and the integrands of X\ 
and Z 2 differ only by an 0(1) multiple, the above analysis also yields that (W^)~ 1 Ii 
is at most of the same order as (H / '^)~ 1 I 2 . Finally, we estimate (W^)' 1 Dk(0) Z 3 , cf. 
(16T20]) - ([6T2T|) . First, we estimate 



exp (±jjf\ dy = e^5-'C>>, 



for some 0(1) constant C' 3 \ Substituting into (I6.2Q[) . then, we obtain 



where 



(^)" 1 J D fc (0)X 3 = £- 1 / 12 
hx(x,z) 




l 3 (x,z) exp 



dA„ 



a 



and n 3 (x, z) = J_(x*) + J-(z) — 2I(x), 



" 3 v ^)F 1 /4( 2 ) 

for some 0(1) constant C 3 . The exponent Il 3 has a sole maximum at the point 
(x**,x**) G d(ir y D2) which is also not a critical point (compare to the maximization 
of Hi in Section [3]). Hence, Theorem Appendix D.l yields 

H 3 (x* , x* 



e 2/3 C 3 exp 



3(x*,x*J exp 



n 3 (x**,x**; 



for some 0(1) constants C 3 and C 3 and where Il 3 (x**,x**) < 2J_(x*). 

In total, then, and to leading order, we obtain the leading order formula 

, Cf d k ( k (x*) /i 2 (0,a;*) 



''OOfc 



for k e 



-1/3 



(6.22) 



Here, we have used that f(x*) 
(cf. (ICTIj) and (1CT20]) ) 



to leading order, while /i 2 is given in ( I3.6P and 



= Ai'(Ai) (Bi^; 



(7, 



2/3 




— and du = - 

irC 3 (N k + A ) 



To derive the desired formula ( 16. ip from ( 16.221) . we note that (cf. fl3.5p -f l3~6lT 

sinh (v / Aq(1 - £*)) 



h 2 (0,x*) 
Hence, (I6.22p becomes 



■0/(0) 



/ A~^ cosh \/Ao 



(6.23) 



1 - v) Cf a 2 /3 f(0) G(x*) sinh (^(1 - x,)) 



l 00k 



a* C 3 (N k + A ) cosh ^A^ 

The desired formula ( 16. ip may now be derived from this equation by recalling ( 12. 5ft and 
the definitions collected in (14. 2p . 
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6. 3. Higher order terms in the asymptotic estimate for a' QOk 

As became evident in the material presented above, certain terms among those we 
estimated are A — dependent, and hence they do not necessarily remain higher order for 
asymptotically large values of A . As we will see in this section, certain terms which 
are higher order for Ao = 0(1) become leading order for Ao S> 1. Apart from that, 
these higher order terms have an important effect even for A = 0(1), as they lead 
to the singularly perturbed problem (I4.38[) for the steady states of the reduced system 

To quantify these terms, we recall from the last section that 

2tt 



(Wg- 1 D k (0) Z 4 + (W^ 1 X 6 + 1 8 = -e 



1/3 _ 



J— ( X J* 



J 4 1 X^ i X * 



By definition of S4, 

^Ay&* 1 X* \, 



1 — '4 \ X* 5 X% 



E ^ 1 — ig (x* X* ) S\ — ig {x* X* j 



(6.24) 



(6.25) 



where S4, Ee, and S§ are expressed in terms of the function hi defined in ( 13. 6ft — see 
(16. 6p . ( I6.8p . and ( I6.10p . respectively. As we saw in the last section, 

E- l 'H 2 D k (U) I (\-v)Cl oT f(0) Ckfa) sinh (yggl - xj) 
^ 4 (J* C 3 v 7 ^ (A fc + A ) cosh v 7 ^' 

At the same time, we calculate 



£ 



/3 S 2 1/2 (1 - u) C 2 C$ al'\ k {x,) sinh (^A~ (1 - x*)) /** /(x) cosh ( v^x) dx 

— J 6 — £ 



1/3X2 
— Is 



2(7, 

(1 - 1/) C? Cf ao /3 a(x*) cosh (y/A^x*) sinh (v^(l - s*)) 



Aq cosh v A 







t o* VAoCoshvAo 

There are two distinguished limits for these expressions, namely, 



- 1/3 ^ 2 Pfc(o) 

£-V3 a 2 



I 4 



eh 



and 



-V^ 2 D fc (0) 



<r* C 3 N k +A ' 

.1/2 (!-»/) (l-g^CgCl^q^) Jo"* f(x)dx 
2a, ' 

. (i-i/)(i-s.)c?c|<rg /3 6fc(x.) 



(l-i/) Cf <rg /3 /(O) Cfc(^) e-v^ 



for A <1, 



(6.26) 



1/3,52 



-l/3 5 2 



^.C 3 V^(JV fc +A ) : 
,1/2 (!-»/) ICfCg^CfcQE.) J_ 

4cr» A ' 

. (i-^)c 2 c 2 CT g /3 c fc (x») 1 

2ct* 



for A > 1, 



(6.27) 



where we have used Theorem Appendix D.2| to estimate the integral appearing in the 
definition (16.81) of S 6 . It immediately follows that eS 8 (x*,x*) <C e 1 ^ 2 S 6 (x*, x*) for all 
A « e' 1 ' 2 . 

Next, we estimate (W^)" 1 !^ + eX 7 in the regime A 1. First, we recall ( 16.1 7ft . 



(W^)- 1 l 5 + £l 7 = £ 



1/3 




E B (y,z) exp 

D 7 V V £ 



dA 
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where E 5 (y,z) = E 5 (y,z) + e 1 / 2 H 7 (?/, z). The functions S 5 and S 7 are expressible in 
terms of the function h\ defined in (13. 5p . see (I6.14p and (I6.16p . respectively. Working as 
for /12 above, we procure the leading order asymptotic relation 

h(x, y) = r f(x) (l - - 9 (x - yj) h 2 (x, y). 

Here, 9{s) = (1 — e _2s )/2 — and hence 8(0) = — while the first term in the right 
member is Ao— independent and hence remains bounded in this regime also. Using this 
expression, we can establish that (W^)~ 1 T§ + Z 7 is at most of order e 4//3 Aq 1 and hence 
higher order. 

Similarly, (I6.18P yields to leading order 

m) X2 ~ £ 5 2F3/4( X ,) 7^' fOT Ao>>1 ' 

where C 2 and C' 2 n are 0(1) constants. Hence, this term is also higher order. The 
term (W^,) _1 Zi can be bounded in a similar way, whereas (H 7 ^,) -1 Dk(0)X 3 is, here also, 
exponentially smaller than all other terms. 

In total, then, and to leading order, we obtain the formula 



a 00k 



' 1 s n2 2/3 _ w s( /(0) sinh(y^(l -X*)) 

V C 3 VAq (iV fc + A ) smh VA 



1/2 sinn ("v/Aq (1 - x *)) jo* f( x ) cos h (\/A~ox) rfx 



2 vAo cosh VAo 

This formula precisely matches (16. 2D . The two associated distinguished limits are 

, (l-^)(l-^)C 2 2 a 2/3 C fc (^) ( /(0) C? ff* /(x) dx ^ 

° 00fc - * l"(iV fc + A )C3 + 2 J' forA o« 1 ' 

and 

a00k ~ * l v "^v^(iV fc + Ao) + ^ A^J' forAo>>1 - 

Note that, in this last formula, the first term in the parentheses dominates the second 
one for all o(l) values of \i (cf. (16. 3p ); the two terms only become commensurate for 
0(1) values of \x. 



7. An asymptotic formula for b' 



mOk 



Finally, we derive the asymptotic formula for b' m0k 

b'mok = -A' k (A ) B, for 0^k,m « s~ l /\ (7.1) 
which has already been reported in (14. ip . We also remark that, here also, this result 



is valid for those values of k for which Cfc(^*) 7^ 0- Theorem Appendix D.l yields an 
(algebraically) higher order result for the remaining values of k. 
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b 'mOk = ^~ 1/6( ^ (1 - v) / f(x)( m (x)u (x)ij k (x)dx 

Jo 

+ eV*6r 1 (l-v) f f(x)Cm{x)u+(x)Ux)dx 
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1 nx 



777- / / f( x ) f{y)Cm{x)Ch{y)^o{x)^hA x )^k+(y) d y dx 

vv i> Jo Jo 

777~ I / f(x) f(y)(m(x)(k(y)uJo(x)ip k: -(x)ip£_(y)dydx 

W</> Jo Jo 



+ 



W f jo 
1 



1 rl 



f(x) f(y) Cm(x) (k(y) wo(x) ip k ,+( x ) $i-(y) dydx 



+e 



f(x) Crn(x) ( k (x) CJ^(x) dx 







(7.2) 



Let Xi, . . . , X 4 denote the integrals in the right member of this formula in the order 
that they appear in it. We will derive the leading order terms in the asymptotic 
expansions of these integrals using Theorem Appendix D.l as in the previous section 
and also for k,m <C e^ 1 ^ 3 . In what follows, we omit the terms Oujq _(1; Xo)wo,+(1; ^o) 
and 9uq _(1; xo)cu f + (l; £o) in (1A.1I) and (1A.2I) . respectively, as one can show that their 
contribution is exponentially small compared to the leading order terms (see also 
Section [3] and EHE]) . 

First, we derive a formula for — {W^)~ 1 D k X2 + (W^) _1 X 3 + 5X4. We write 

X, = 



f(x) Cm(x) u (x) ^k,-( x ) dx f(y) Ck(y) %Jy) d v 



Jo 



f(y)Ck(y) ^Jy) d y, 



where we have used that ip k - — C^ 1 a 1 ^ 3 loq in a neighborhood of the origin, that 
u is exponentially small outside this neighborhood, the identity 1 1 c*^ 1 1 2 = 1; an d f !2.5j) . 
Employing (IC.12p . next, we obtain 



.7/12 



/(Q) C 2 f 1 f(y) Ck(y) 

^C.o-T Jo F^{y) 



cxp 



v 7 ^ 



dj/. 



Substituting for D k from ( 1C 19|) . we obtain 



(W^)- 1 D k l 2 = e 



-1/12 



;(?/) ex P 



My) 



where we have defined the functions 

„ C 2 f(0)d k f(y)( k (y) 

-2{y) - 



and n 2 (y) = J-{y). 



(7.3) 



(7.4) 



Next, we change the order in which integration is carried out in X3 and use ( lA.ljl and 
(IC.12p -f lC.13P to rewrite this integral as 

(W^-'Xs = {W^Y 1 jf ^ f(x)( m (x)u;o(x)^ + (x)dx^ f(y)( k (y)i>+_(y)dy 
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.5/12 



£3(2/) exp 



My) 



dy, 



where n 3 (y) = II 2 (y) and 



My) 



c.c^l" ( f y f(x)C m (x) J \ f(y)Ck(y) 



4vr 3 / 2 W, 



VW) 



dx 



Finally, using ( 1A.2|) and renaming the integration variable x into y, we obtain 



eTa — e 



-13/12 



My) ex P 



n 4 (y) 



dy, 
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(7.5) 

(7.6) 

(7.7) 



where 



My) 



CiC 2( x 1/3 /(y)Uy)G(y) 



2^ FV'(y) 
Combining f l7.3p - fl7.8p . we obtain 

-(^)- 1 J D fc X 2 +(^)- 1 X 3 +X 4 = -e 



and n 4 (y) = My) = My)- (7.8) 
Mv)' 



-1/12 



H 2 (y) exp 



v 7 ^ 



where, to leading order and uniformly over [0,1], 

C 2 d k (k(y)f(o)f(y) 



My) = My) 



Regarding X 1; we use ( 1A.2|) and flC12[) (lC13|) . to write it in the form 

ni(x,y)' 



(W^)- 1 l 1 =e"/ 12 J j D M^y) exp 



dy, (7.9) 



(7.10) 



(7.11) 



(7.12) 



where D = {(x, y)|0 < x < landO <y< x}, Ui(x,y) = J+(y) — 2I(x), and 

„ ( C^al" f(x)f(y)Ux)Ck(y) 

^ 1{,V) 4vr 3 / 2 ^ y/F{x)F^{y) ' 

First, we estimate -{W^)~ l D k X 2 + (W^)^ 1 X 3 + el 4 , cf. (I7l^ - (I7TT0"]) . The 
exponent n 2 assumes its maximum at the interior critical point G (0, 1), and hence 
Theorem Appendix D.l yields 

-(^)- 1 J D fc X 2 + (^)" 1 X 3 + eX 4 = -e 1/A 

Here, 



'2tt 



(e- 1 l 12 5- l M^))=-e 1/& 5- 1 C 2 . 



V2C 2 o-l /3 f(0)f(x*)( k (x«) 



C 9 



dC 3( rl /2 (N k + A ) 

Next, we estimate Xi, cf. (I7.lip - fl7.12p . The exponent TEi assumes its maximum at the 
point (x*,x*) G 3D which is not a critical point of TIi (compare to the maximization of 
n 4 in Section [3]). As a result, Theorem Appendix D.l| yields 

my 1 !, =e^C[ (e*' vl 5- 1 E 1 (z.,z.))=e 7 '*8- 1 CZ 

to leading order, and with C[ and C" being 0(1) constants. 
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for m, k <C e 



1/3 



Formula (17. ip now immediately follows. 
8. Discussion 

As argued in the Introduction, there are two contextual themes central to this article. 
The first one relates to understanding the nonlinear, long-term dynamics of small 
patterns of DCM type generated through the linear destabilization mechanism identified 
in [25] . The second theme concerns the development of a concrete approach to studying 
the dynamics generated by the (rescaled) PDE model (jl.5p near a linear destabilization 
but beyond the region of applicability of the center manifold reduction. In this article, 
we have reported significant results (outlined in the Introduction) touching on both 
themes. These results, in turn, inspire further investigation within this dual context. 

Regarding our first focal point, and in view of our discovery that the bifurcating, 
small-amplitude, DCM pattern undergoes a Hopf bifurcation, the central question is 
naturally what happens beyond this secondary bifurcation. This question can be 
answered by the methods developed here, as it is in principle possible to deduce 
analytically the sub- or supercriticality of the Hopf bifurcation undergone by (14. 8ft . 
The numerical simulations of [15] indicate that this bifurcation may be only the first 
of a cascade of subsequent period-doubling bifurcations leading to a region of spatio- 
temporal chaotic dynamics and throughout which the phytoplankton profile maintains 
a DCM-like structure. There is, of course, no a priori reason for this cascade to occur 
entirely within the regime A — A* = 0{e) covered by our analysis here. In fact, the 
simulations of [15] suggest that, for the parameter combinations considered there, this 
is indeed not the case. On the other hand, our analysis is able to determine regions 
in parameter space where this cascade can or cannot occur (for instance, in the event 
that the Hopf bifurcation turns out to be subcritical). Moreover, the possibility that 
there exist regions in parameter space where the entire cascade is within the reach of 
our asymptotic methods cannot be excluded. A similar question concerns, naturally, 
the origins and fate of the second DCM pattern identified in Section 14.51 

These last remarks bring us to the second theme. The approach we developed 
here will be used — and if necessary extended — in forthcoming work investigating 
the remaining issues pertaining to our linear destabilization results in [25] — namely, 
determining the nonlinear behavior associated with the destabilization of BL-type. Our 
analysis in [25] strongly suggests that, for realistic choices of the parameters pertinent to 
shallower water columns {e.g., estuaries and lakes), patterns of benthic layer (BL) type 
are equally relevant to the dynamics generated by (II. ip as the DCM patterns considered 
here. In fact, preliminary numerical simulations strongly suggest that co-dimension 
two-type patterns combining DCM and BL characteristics play an important role in 
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the region where the trivial state is unstable. From a mathematical point of view, the 
co-dimension two point may also be seen as an 'organizing center' for the more complex 
behavior exhibited by the system studied numerically in [15]. That is, the cascade of 
period doubling bifurcations reported in [15] may be based on the presence of that co- 
dimension two point. In view of that, the derivation and analysis of an extended reduced 
system for parameter values valid within an 0(e) neighborhood of that point may prove 
highly engaging. 

The same methodology can also be applied to extended models. A natural extension 
of (11. ip is a multi-species model, i.e., a model similar to (11. ip in which several 
phytoplankton species compete for the same nutrient. At the linear level, the species 
evolution decouples [25]. Nonlinear coupling, however, is present through shadowing 
(light limitation) and nutrient uptake (nutrient limitation), and hence the presence of 
every extra species affects the life cycle of each species. Reaction-diffusion models of 
this sort for eutrophic environments — i.e., in the presence of an ample nutrient supply — 
have been developed and investigated both numerically [2] and theoretically [8]. The 
oligotrophic case, on the other hand — where these multi-species models are coupled to 
a PDE for the nutrient — has so far only been investigated numerically [T5] . 

Another natural, if not outright necessary, extension is the inclusion of horizontal 
spatial directions. Plainly, the dynamics generated by (II. ip will be strongly influenced 
by the flow in directions perpendicular to the one-dimensional water column considered 
here: oceanic currents are bound to mix neighboring water columns and thus also enrich 
the collection of emerging planktonic patterns. Finally, and as already described in the 
Introduction, we are currently studying the simplified model problem (11.211) through 
which we hope to understand the applicability and limitations of the general method 
developed here. This approach may also serve as a first step towards obtaining a rigorous 
validation of our method. 
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Appendix A. An asymptotic formula for ujq 

The formula for the principal part in the asymptotic expansion of Uq reads 



cf. [25], where x , C\, C 2 , F, o~ and 9 have been defined in (11.151) . (J27TTJ) , (13. 3p . and 
(I3.12p . We remark that C\ is a normalizing constant ensuring that ||wo|| 2 = 1- (This 



u (x) 




(A.l) 



u ,±(x;x ) = exp ± 



Wo±(x; x ) = E(x) uj 0: ±(x; x ) = exp 
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factor does not appear in the formula for uq we give in [23] , since Uq was not normalized 
there.) Also, 

where / has been defined in (I2.16p . An asymptotic formula for Uq = Euq is readily 
derived using (1A.1I) above, 

( e~ 1 / 6 a 1 /e C l eV^-M(A 1 (l-x^x)), for x G [0,x ), ^ 

0Jq{x) ~ < e -i/i2 C c CT i/3 (A.2J 
6 2^f^/?{x) [ u o^( x ' x o) + 9 x ) u^ + (x; x )] , for x G (xq, 1], 

where we have defined the functions 

' J±{x) 

with J± as in (I2.16p . We remark that J_ becomes maximum at the well-defined point 
x* G (0, 1) — the location of the DCM, see ( 12.181) — whereas J + increases monotonically. 
Also, the terms involving co>o,+ in (1A.1I) and 0Jq + in (1A.2I) are exponentially smaller than 
the terms Uq (x) and u^_(x), respectively, everywhere except for an 0(y/e)— region of 
x = 1. Indeed, for all x < 1, 

J+(x) - 2/(1) = J_(x) - 2(1(1) - I(x)) < J_(x). (A.3) 

In particular, Hu^Hco can be bounded by an 0(e~ 1 / 12 5~ 1 ) constant, where 6 = 
exp(— J_(x*))/\/e is an exponentially small parameter (cf. (I2.15jl ). 

Appendix B. An asymptotic formula for rj Q 

We recall that rj is the solution to the boundary value problem ( 12. 6p . 

£d xx r} - \ rj = -ef^fufi, where d x r} (0) = r] (l) = 0. 
Recalling that Ao = sAq in our bifurcation analysis, we find that 

d xx 7] - A ?7o = -i' 1 /^, where d x 7] {0) = r] (l) = 0. (B.l) 

The functions T)q ±(x) = e ±v ^ x form a pair of fundamental solutions to the homogeneous 
problem. Using variation of constants, then, we obtain a special solution to the 
inhomogeneous ODE, 

Vo,s P (x) = (2£ v / Ao)~ 1 [To (vo,+f ; x) Vo-(x) - T (rj -f ; x) Vo,+( x )} ■ 
Here, we have defined the family of functionals 

T n (■ ; x) = / -(s) UJn( s ) ds, parameterized by x G [0, 1] and n > 0. (B.2) 
Jo 

The solution to (IB.ip is, then, 

Vo (x) = [C+-(2£ v ^)- 1 r (r /0i _/;x)]r /0 , + (x) 

+ [^ + (2£ v ^)^r (r /0i+ /;x)]r /0i „(x). 1 ' ' 
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'A - VAo C~ = 0, 



Imposing the boundary conditions for r] and using the identity T Q (• ; 0) = 0, we find 
that the constants C~ and satisfy the linear system 



2VA C+-r (%,_/; 1) 
the solution to which is = C~ = C^/(2^\/Ao), with 

r (%,-/ ; i) - r (vo,+f ; i) vo-i 1 



2£^%C- + T (770,+/ ; 1) 



2 cosh yAo 
Thus, (1B.5|) becomes 

^o(x) = (2£ V / A^)~ 1 ^cosh^v^Ao^) + r o ; *0 fto,-0*0 - T (%,_/; a) t^-hO*)] .(B.4) 

Further employing the definition (IB. 21) . we calculate 
To (Vo,+f; x)t]q_{x) - T (rjo-f]x)rj , + (x) 

[vo,-( x )vo,+(y) - vo,+( x )vo,-(y)\ f(y) ^oiv) d v 



sinh ^A (x-y)j f(y) u£(y) dy 
2T (sinh(y/A (x--))f;x] 



Additionally, 



r (vo-f ; i) ??o,+(i) - r (770,+/ ; 1) vo-(i) 



2 cosh \/A 







1 



2 cosh yAo jo 
1 



cosh vA 



Jo 



[VoAv) - %,+(!/) %,-(!)] /(?/) w£(y) dy 

sinh (v / A "(l - y)j f(y) u£(y) dy 



cosh v A 
and hence (1B.4[) becomes 

r) Q (x) ~- 



^-=r (sinh (VXo"(l -■))/; 1 



cosh (VAflx) / / / — , \ 
_i^lr (sinh (^(1 -))/;! 

-r (sinh (\ZA^(sc -•))/; 3 

COsh (-v/AnX) Z" 1 / / , \ I 

COS h ^ y sinh " y) ) M<4to)*u 

sinh ( y/A^(x - y) ) /(j/) dy 



(B.5) 



To estimate ||^o||oo over [0, 1], we first show that 770 is positive and that it assumes 
its maximum in an 0(e 1 ^ 4: ) neighborhood of x*. First, an estimate based on (1B.5[) 
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f 1 cosh (x/Aq'x) ( i — / x\ ( i — , 

/ ^ ' sinh V A o(l~ y) -sinh (x/A (x-y) 

Jo coshVAo v J v 

sinh (v / Aq(1 — x)) 



dy 



cosh V A y) f(y) u+(y) dy > 0, 



£vA coshvAo jo 

for x G (0,1). To locate the maximum, we differentiate both members of ( 1B.5P and 
obtain 



£d x rj Q (x) 



sinh (a/Aox) 
cosh y/Ko 

— / cosh 



sinh v A (l - y) f(y) Uo(y) dy 



Ao(^-y)) f{y)uJo(y)dy 
sinh (y/7^(x - y)) 



(B.6) 



Theorem Appendix D.l can be used to yield the principal part of the two integrals in 
this formula, whereas the term proportional to ojq can be estimated via (1A.2j) . For the 
values of A we are interested in, the localized term in either integrand is Uq, while the 
A — dependent terms vary on an asymptotically larger length scale. Thus, 

sinh(VA^x) f l 

o x r]o{x) = -=J- I 

£coshvA Jo 

to leading order and for x < x* and \x — x*| ^> e 1 / 4 , since the second and third terms in 
the right member of (IB. 61) are exponentially small compared to the first one. Similarly, 



sinh 



Ao(l-y) f(y)uj+(y)dy>0, 



d x r]o{x) 



cosh vAo 



sinh ( v A x 



— cosh \J A 



sinh 
cosh 



'Ao(l-y)J f(y)u+(y)dy 
/Ao(z-y)) f(y)uo(y)dy 



for x > x* and \x — x*\ 3> £ 1//4 , since the second and third terms in the same formula are 
of the same asymptotic order and the third one is exponentially smaller. Changing the 
upper limit of the second integral to one (and thus only introducing an exponentially 
small error) and combining the two integrals, we find 

cosh (v / A^(l - x)) rl 



d x rj (x) 



cosIi^a/A^?/) f{y)uo(y)dy < 0, 



i cosh a/Ao 

Since r] G C 2 (0, 1), now, it follows that r]' (xi) = at a point X\ such that \x* — X\\ 
0(s 1//4 ), as desired. Hence, we can now use ( 1B.5I) to estimate further 



cosh (\/A^xi) 
t \f~Ko cosh \f~Ko Jo 



sinh ( vAo(l - y) ) f(y) Uo(y) dy 



Using our asymptotic estimate on X\ and Theorem Appendix D.l , we find 

1/6 ! cosh (v/A^x*) sinh ( V / A^(l - x*)) 



1*701 



/ Aq cosh a/Aq 
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for some A — independent, 0(1) constant C. Since the A — dependent quantity in the 
bound above remains bounded by an 0(1) constant also for A ^> 1, we finally conclude 
that 1 1 ?7o 1 1 oo can be bounded by an 0(e 1//6 <5 -1 ) constant. 

Appendix C. Asymptotic formulas for ip n , n > 

The function ip n is the solution to the boundary value problem 

£ d xx ip n + {f{x) - £ - v - v n ) rj) n = -e£~ l fE ( n , where Q (if> n ; 0) = Q {ip n ; 1) = 0, 



cf. ( 12. lOD . Here, Q (ifj n ; x) = ip n (x) — ^/ e / v d x ip n (x) and we recall that 

C n {x) = V / 2cOs( V / iVnX), (C.l) 

see (I2.4p . Recalling also the definitions F(x) = /(0) — f(x) and A* = /(0) — £ — v, as 
well as that A* = A + s^Vo by (11.131) . we write 

/( x )_^_ v = Ao-F(x)+£ 1/; Vo, 

with yU = o"q/ 3 I Ai| + 0(e 1 / 6 ). Finally, since Ao = eA and v n = —eN n , we may rewrite 
1 12. 10p in the final form 

ed xx ^ n - [F{x) - e 1/3 /i -e(N n + A )] ^ n = - £ ^ /Cn , (C.2) 

together with the boundary conditions (-0 n ; 0) = Q {ij) n ; 1) = 0. In what follows, we 
derive asymptotic formulas for ij) n and for values of n satisfying n <C e^ 1 ^ 3 . In that case, 
e(N n + A ) -C e 1 ^ 3 — recall our assumption that A <C e~ 2 l' 3 in Section 12. 3 .21 — and hence 
this term is perturbative to e 1 ^ 3 ^. Hence, we may write 

e 1/3 fi + e(N n + A ) = F(x n ), where x n = x (l + o(l)) (C.3) 

is a turning point for ( 1C.2I) . Then, (1C2|) becomes 

e d xx *p n - [F(x) - F(x n )} i, n = -ei-'fE ( n , (C.4) 

equipped with the boundary conditions ( 12. 1011 . The solution to this boundary- value 
problem may be found by variation of constants, 

i/j n (x) = [C+ - (iW^G^x)] ijj n , + (x) + [C£ + (IW^G+ix)] il>n,-(x). (C.5) 

Here, ip n ,± is any pair of fundamental solutions to e d xx ij; n = [F(x) — F(x n )] ip n and 
= ipn-dxipn^ — ^n,+d x 4> n - is the associated Wronskian. (To derive the result 
above, one needs to show that is constant. This is plain to show by using the 
identity d x W^(x) = 0, for all x G [0, 1], which follows from the definition of and the 
ODE that ip± satisfy.) Further, 



G±(x)= / f{y)Uy)<±{y)dy, (c.6) 

JO 

where ipn± = Fip n ,±- Using ( 1C.5I1 . we further obtain 

dMx) = [CJ - (iW^G^x)} d x ^ + {x) + [C£ + {£W^)- l G + {x)] 9^ n ,_(x), 
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and thus the boundary conditions yield the system 



■GJ1) 



The solution to this system is 

1 



Q (Vw ; i) + 



c: 



-G + (l) 



g (V„,_ ; i) = o. 



where 



r7+ 



A/, £ C0n - 5 0) and C, 



1 



G_(l)0(^ + ;l)-G + (l)0(V> ni _;l 



£ ; 0) g Wn,- ;i)-g ; o) Q (i> n ,+ ; l) 

Thus, also, (1C5|) becomes 



D^g{^ n>+ ;0), (C.7) 



(C.i 



^(s) = (W,/,)- 1 [r_(x)^ n ,_(x) -r+(x)V„,+ (a:)], 



(C.9) 



where 



r_(x) = G+(z) + D^g {il> n>+ ; 0) (CIO) 
r+(x) = G_(x)+Z^£ (</>„,_ ;0). (C.ll) 

These formulas hold for an arbitrary pair ^/> n) ± of fundamental solutions. Working 
as in [25] , where the problem was considered in detail in the absence of the perturbative 
term e(N n + A ), we can derive the following leading order formulas for a specific pair 
of solutions ipn,±'- 



1pn,-(x) 



e^a'^Ai (Aril 



x, 



-1/4 



C 2 



1/6 ff -V6 



-1/4 



Bi (A 1 {l-x^ 1 x)) 



for x G [0, Xo), 
for x G (x , 1], 

for x G [0, x ), 
for x G (xq, 1]. 



(C.12) 



(C.13) 



Here, we have used that x n 
earlier, leads to 



xo + o(y/e). The identity d x W$ = 0, which was reported 

-Ai\A 1 )Bi(A 1 ) = lim W4x) = 1/tt > 0, (C.14) 

for all x G [0, 1] and for this particular pair. (To calculate the limit, we used the 
asymptotic expansions of Ai(%) and Bi(x) as x ~ ► 00 — see ; e -9-, PJ-) Next, we simplify 
the formula (1C.8|) by investigating the asymptotic magnitude of the terms in its right 
member. By definition (12.101) . 

g (Vy± ; o) = ^,±(o) - v / ^(^^n,±)(o). 

Equations (JC^D and (ICT2|) - dCT3|) yield 

g(^n,-;0) = -£ 5 / 6 a - 5/6 Ai / (A 1 )(iV„ + Ao) + 0( £ 7 / 6 ), 
£ (^n,+ ; 0) = e l/Q a - 1/6 Bi(Ai) + O^ 3 ). 
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(Here, we have Taylor expanded Ai(yl 1 (l — x$ x)) around its zero x = 0.) Next, 



1/4 (1 =F y/<Ti/v) C± ( , 1(1) 



g& n ,±;i) -g 1/4 V \;l ' " exp^±^j, (c.15) 

recall A3. 12[) . These formulas imply that Q (ip n ,+ '■, 0) Q (Vv,- j 1) is exponentially smaller 
than Q (if) n ,~ > 0) Q (i } n,+ ] 1), and thus 

^,(1) ^(1) - G-(l) ^ g(^;l) 

* ff(A,-;o) 5 (*,+ ;i) 

and down to exponentially small terms. Next, the relative asymptotic magnitudes of 
the terms in G_(l) — _D n (l)G + (l) may be derived using the definitions (12.101) and 0C6p 
together with Laplace's approximation (cf. Theorem Appendix D.l ). One finds that 
G_(l) is dominated by exp(e~ 1 / 2 J_(x*)), whereas D n (l)G + (l) by exp(e~ 1 / 2 J_(l)), and 
hence the latter is exponentially smaller than the former. Hence, 

Q Wn,~ ; 0) 

It follows, then, that 

r_(x) = G + (x)-D n (0)G_(l) and r+(x) = G_(x) - G_(l), (C.18) 
and down to exponentially small terms. Here, 

y (V>n - ; 0) 

where (recall (1C.14I) ) 

Combining this formula with ( 1C.9|) . we find 

^ n {x) = (iW^)- 1 [G + (x)ij> n ,-(x) - G_(x)^ n!+ (x) + G_(l) & n ,+(x) - D n (0)^ ni _(x))} 
= (IW^y 1 [(G + (x)-D n (0)G_(l))V„,-(x) + (G_(l)-G_(z))V„ J+ (x)] 



4>n,-(?)[ f(y)Uy)K,+(y) d y- D n(o) f(y)Uy)4>Z-(y)dy 



/ f(y)(n(y) i>nAy) d y 



(C21) 



Appendix D. Asymptotic approximation of integrals 

Appendix D.l. Localized integrals 

Our main tool in this section will be Laplace's method and, in particular, the following 
three theorems based on [2U Ch. II, VIII, IX]. 
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Theorem Appendix D.l Theorem IX. 3]) Letn 6N,Dc R" be a domain with 

piecewise smooth boundary dD, and Uq G D. Let, also, the functions II G C 2 (D, R) and 
H G C(Z), R) satisfy the conditions 

(a) inf II(u) > II(u ), for all 5 > 0, 

D-B(uo;5) 

ft) a (£> 2 n( Mo )) C R+, 

(c) the integral Xd{\) := J' ''J E(u) e~ xn ^du converges absolutely 
for all sufficiently large A. 
(Here, D 2 U denotes the Hessian matrix of II.) Then, 

oo 

T D {\) ~ e- An « Yl Ck A" (fc+n/2) (A — > OO), 

where one may derive explicit formulas for the constants {ck}k- In particular, 
/2vr\ n/2 - (un) e - An ( u °) 

(I) lD(X) ~(t) ^ CTM ' !/ -^W- 

(II) Z D (A) ~ (-J C e- An ^, if u Q eD and E(u Q ) = 0, 

(III) T D {\) ~ — - ° J n?TT/ . % G <9D, S(«o) + 0, and Dn(« ) = 0, 

\ A / 2y det -D^I^Mo) 

/2vr\ (n+1)/2 » (unl e - An(uo) 

(IV) Z D (A)~(^ T J "l^etJ ' ^ U ° e HM ^ °' ^ DU{U0) ^ °' 
as A — >■ oo, /or some constant Cq which is at most 0(1) with respect to A and under 
the assumption that dD is smooth around u in the cases where Uq G dD. Here, J is a 
matrix related to D 2 U(uo) and to the local characteristics of dD around uq. 

Theorem Appendix D.2 Let a < b and uq G [a, 6]. Let, also, the functions 
II G C 2 ([a,6],R) and H G C([a, 6],R) satisfy the conditions 

(a) inf II(w) > n(« ), for all 5 > 0, 

[a,b] — B(uo;S) 

(b) the integral X(A) := / S(u) e _An( - w W converges absolutely 

J a 

for all sufficiently large A. 

Then, 



oo 



X(A) ~ e" An ^ ^ c k A" fc / 2 (A ^ oo), 



fe=i 



where one may derive explicit formulas for the constants {ck}k- In particular, as A — > oo, 
(I) X(A) ~^T7^^==, n G(a,6) and E(u Q ) ? 0, 

A 1 / 2 ^/n "(u ) 
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e -An ( n ) (~> ) - E ' ( tZ) U0) ) 

(II) 1(A) — J -, if u e(a,b) and E(u ) = 0, 

v > ^ > A 3/2 v^[n"(n )] 3 / 2 

(HI) i(A) — yV-^U, if u e{a,b}, ~(u )^0, and n'(u ) ^ 0, 

a |n'(uo)| 

-An(« ) .fj?~( Un \ 

(IV) 1(A) ~ -^75- -^===> % G {a, 6}, S(«o) ^ 0, and il'(n ) = 0, 

(V) 1(A) ~ : -^^t^ A ^ ) «oH ,S(n ) = 0, and U'(u ) ? 0, 

(VI) 2(A) tf ,S(«o) = 0, and #(«„) = 0. 

Theorem Appendix D.3 Let D C R 2 &e a two-dimensional domain with piecewise 
smooth boundary 3D and uq £ 3D. Let, aiso, i/ie functions U £ C 2 (Z),R) and 
H £ C(l), R) satisfy the conditions 

(a) _ inf II(m) > n(w ), /or a// 5 > 0, 

D-S(m ;<5) 



A 


|n'(«o)r 


e -An(«o) 


0rH(it o ) 


A 1 / 2 


v/2n"(«o) 


e -AII(jio) 




A 2 


[n'K)F 


e -An(uo) 


±S'M 


A 


n"(«o) ' 



(b) the integral Id (X) := J " 'J ^( u ) e Xn ^du converges absolutely 

for all sufficiently large A. 

Assume, further, that 3D has a corner at Uq and, in particular, that dD is given (locally 
around Uq) by the curves k(x,y) = and h(x,y) = with Dkiuo) x Dh(uo) ^ 0. Let 
the vectors v k and Vh satisfy 

Vk -L Dk(u ), v h JL Dh(u ), and \\vk x %|| = 1. 

and Vh can be selected to further satisfy the conditions 

n fc := (v k , Dn(«o)) > and II h := DU(u Q )) > 0, (D.l) 

i/ien 



oo 



fc=0 



where one may derive explicit formulas for the constants {c k }k- In particular, 
(I) 1 D (X) ~— 1 0) - if E(u o )^0, 

u m ] a 2 2n fc n, v ^+n 2 1 

an r m 1 (n fc S fc + n h 5 fe ) e~ An ("°) 

II In A ~— . , C(nn) = 0, 

as A — > oo. .Here, E k := (v k , DE(u )) and E h := (vh, DE(u )) , compare to W.l\) . 
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Appendix D.2. Oscillatory integrals 

Theorem Appendix D.4 Let a < b, E G C([a, 6],R), and $ G C 2 ([a, 6], R). Assume 
that 

= $( a ) + (t - a) and $'(£) > 0, /or a// 1 G [a, 6] and with $i(a) 7^ 0. 

Then, the integral X(A) := J a b S(t) e lA *^ /ias the following asymptotic expansion: 

1(A) ~ [^^(O) eiA * (a) " ^ (fc) ($(6) - $(a)) e iA * (&) ] lj) (A ->■ 00), 

where we have defined the function 

Mr) =S(f (t)K(t). 

i/ere, r(£) = $(£) - $(a) or, equivalently, t(r) = $ _1 ($(a) + r). 
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